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Abstract 


Improving Particle Confinement in Inertial 
Electrostatic Fusion for Spacecraft Power and 
Propulsion 


By 

Carl C. Dietrich 

Fusion energy is attractive for use in future spacecraft because of improved fuel energy 
density and reduced radioactivity compared with fission power. Unfortunately, the most 
promising means of generating fusion power on the ground (Tokamak based reactors like 
ITER and inertial confinement reactors like NIF) require very large and heavy structures 
for power supplies and magnets, in the case of magnetic confinement, or capacitors and 
lasers in the case of inertial confinement. The mass of these reactors and support 
equipment is sufficiently large that no existing or planned heavy-lift vehicle could launch 
such a reactor, thereby necessitating in-space construction which would substantially 
increase the cost of the endeavor. The scaling of Inertial Electrostatic Confinement (IEC) 
is such that high power densities might be achievable in small, light-weight reactors, 
potentially enabling more rapid, lower cost development of fusion power and propulsion 
systems for space applications. 

The primary focus of the research into improving particle and energy confinement in IEC 
systems is based on the idea of electrostatic ion focusing in a spherically symmetric 
gridded IEC system. Improved ion confinement in this system is achieved by the 
insertion of multiple concentric grids with appropriately tailored potentials to focus ion 
beams away from the grid wires. In order to reduce the occurrence of charge exchange 
and streaming electron power losses, the system is run at high vacuum. This 
modification to the usual approach was conceived of by Dr. Ray Sedwick and 
computational modeling has been conducted by Tom McGuire using a variety of custom 
and commercial codes. 

In this thesis, a semi-analytic model of the potential structure around a multi-grid IEC 
device is developed. A 1-D paraxial ray ion beam envelope approximation is then used 
along an equatorial beamline and the assumed beam density is gradually increased until 
an effective beam space charge limit is reached at which point the potential fusion output 
is calculated. Significant use of the commercial particle-in-cell code OOPIC was made, 
and its ability to predict multi-grid IEC confinement properties is evaluated. An 
experiment was built to confirm the effectiveness of the multiple-grid structure to 
improve ion confinement times. It is shown that the multi-grid IEC can improve ion 
confinement time over the conventional, 2-grid IEC device. The PIC predicted ion 
bunching mode is also seen in experiment. 
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1. Introduction 


1.1. Background 

Inertial electrostatic confinement (IEC) fusion is often credited as being conceived of by 
Philo T. Farnsworth, the inventor of television [20]. In IEC devices, fusion ions are 
electrostatically accelerated through a (nearly) spherically symmetric central focus point. 
At the focus point they are at sufficiently high energy to overcome their mutual 
electrostatic repulsion that some small fraction of the ions will fuse and release energy in 
the form of high energy fusion products. Figure 1 illustrates the practical simplicity of 
IEC fusion. 



Figure 1: Schematic drawing of conventional, 2-grid IEC device 
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The earliest published work on the concept was done by Salisbury in 1949 [1], then 
expanded by Elmore, Tuck and Watson in the late 1950s[2, 3]. There has been sporadic 
work on IEC fusion for the past 50 years [4-350]. Experimental programs are relatively 
cheap to develop, and they have been shown to be a detectable neutron source given 
sufficient input power. The largest fusion output to date was achieved by Hirsch in 1967 
with an output of nearly 10 10 neutrons/second with Deuterium and Tritium [18]. It is 
widely accepted that it is not difficult to make fusion reactions in an IEC reactor, but the 
best efficiency to date as measured by the fusion power out divided by the electrical 
power in (Q) is approximately 10' 5 . Because of its poor efficiency, IEC has never been a 
priority in the quest for a fusion energy source. 

However, IEC is appealing in its simplicity, and it has consequently gained a substantial 
amateur following from people who enjoy the idea of having a fusion reactor in their 
garage. 



Figure 2: IEC "Fusor" built by high-school student Brian McDermott (photo credit: Brian 
McDermott, reproduced with permission) 
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These devices have developed a strong following in online forums and bulletin boards. 
Figure 2 is a photograph of an operating IEC fusion reactor built by high school student, 
Brian McDermott. Many other amateur scientists have built similar “fusor reactors” 
[361], 

Theoretical systems studies have identified significant barriers to the practical 
implementation of IEC fusion systems for power production [175], but significant interest 
remains because of the simplicity, the attractive scaling, and the low mass of the IEC 
fusion reactor system. This low mass is particularly of interest for space-based power 
and propulsion systems. While there are designs for many types of magnetically 
confined and inertial fusion reactors that are much closer to achieving net power output, 
these concepts are sufficiently massive that tremendous monetary resources would need 
to be allocated for development of that type of reactor in space. The scale of such an 
undertaking would likely surpass the scale of the International Space Station (ISS) - 
requiring multiple launches and in-space assembly. It is therefore valuable to work on 
the challenges associated with making IEC fusion more efficient because of the potential 
to put an entire reactor in space with a single launch. The potential for a more practical, 
small-scale, implementation of fusion power for spacecraft systems is the fundamental 
driver for this research. 

The primary limiter of the efficiency of IEC fusion systems is the extremely short energy 
confinement time, the embodiment of which is the short charged particle confinement 
time in the system. Under typical IEC operating parameters, electrons stream out from 
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the central cathode and are collected on the outer anode, and ions typically make only 10 
passes through the cathode before being lost by either impacting one of the grid wires or 
undergoing a charge-exchange reaction with background gas in the system [357]. In a 
typical system, the probability that any given ion will be lost to a process other than 
fusion is approximately 100,000 times more likely than losing an ion to a fusion reaction. 

Previous work [357] identified the need to operate an IEC fusion device in a regime 
where the ion loss probability is within the fusion energy gain to input power ratio of ion 
fusion probability (a factor of order 10) in order to approach Q=l. The improvement of 
charged particle confinement in IEC systems is therefore of paramount interest. Sedwick 
and McGuire proposed a method of improving ion confinement times by electrostatically 
focusing ion beams to keep them away from cathode grid wires. McGuire also identified 
the need to operate the system in the ultra-high-vacuum (UHV) regime in order to avoid 
charge-exchange and collisional scattering losses with background gas. With an 
appropriately hard vacuum and a properly focused recirculating ion beam, typical ion 
confinement times were computationally shown to be improvable by three orders of 
magnitude - from 10 passes to 10,000 passes dependent on the background pressure. 

This improved ion confinement has the potential to yield dramatic improvements in 
system efficiency [357]. 

The practical implementation of this type of improved confinement requires the 
development of an IEC reactor with 3 or more independently biased, concentric spherical 
grids (multi-grid). The potentials on these grids would be set so as to provide a confining 
electrostatic field. By minimizing the rate of ion-grid impact, the electron current will 
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also be reduced because the majority of streaming electrons are emitted as secondary 
electrons when ions collide with the cathode. Electrostatic ion confinement thereby has 
the potential to reduce the two largest energy sinks in IEC systems. 

While focusing ion beams is nothing new to physicists, and there are many tools to 
predict ion beam behavior in the presence of focusing fields, the development of new 
models was required in order to predict the behavior of ions in this type of recirculating 
ion trap. The development of some of those new models and the construction and 
experimental evaluation of the first multi-grid IEC confinement experiment is the subject 
of this thesis. 

1.2. Rationale Behind Current Investigations 

Previous work computationally identified the potential to improve ion confinement in 
IEC devices by implementing a multiple-grid configuration, but hardware based tests 
were needed to experimentally validate the predicted improvement. 

1.2.1. Ion Confinement 

The goal of this research is to compare a conventional, 2-grid IEC device with a 5-grid 
IEC device to evaluate the potential of the multi-grid approach to improve ion 
confinement. Based on the computational modeling, significant increases in ion 
confinement were expected with properly tailored potentials in the 5-grid device. The 
major ion loss mechanisms in a conventional, 2-grid IEC device are chaotic ion 
trajectories and background gas interactions. Operating the 5 -grid device at very low 
background gas pressure and with a focusing field structure is expected to improve ion 
confinement. 
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I.2.I.I. Unconfined ion trajectories 


Particle-In-Cell (PIC) simulations conducted by McGuire and repeated with modification 
by the author indicated that the typical 2-Grid IEC device would be subject to chaotic ion 
trajectories due to both field asymmetries from the presence of a feed through and the 
highly defocusing nature of the vacuum electrostatic field structure surrounding the 
cathode. This defocusing field curvature in the region of the cathode grid wires could 
cause a recirculating ion beam to either disperse resulting in spatially chaotic ion 
trajectories or impact a cathode grid wire after only 2-3 passes through the device core. 



Figure 3: OOPIC model of 2-grid ion trajectories with field asymmetry due to the feed-through 
Figure 3 shows a particle in cell model of ions in a 2-grid IEC device with a feed-through 
field asymmetry. When multiple, concentric, independently biased wire grids are used as 
shown in figure 4, the vacuum electrostatic field structure was shown to be capable of 
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confining the ions away from the cathode grid wires for many more passes through the 
core than in the conventional, 2-grid design simulated above. With the confining field 
structure of the multi-grid configuration, significant improvements in ion confinement 
were shown to be possible even with a small, off-axis ion injector that would act as an ion 
sink as well as a source. 



Figure 4: Multi-grid with asymmetric feed-through and off-axis, absorbing ion injector 

These early studies strongly indicated that uncontrolled ion trajectories in conventional, 
2-grid IEC devices could be a major ion (and energy) loss mechanism and supported the 
hypothesis that the use of multiple grids could significantly improve ion (and energy) 
confinement times in IEC devices. 
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1.2. 1.2. Scattering and charge exchange with background 


In the regime where most conventional IEC devices are operated for maximum measured 
reaction rates, uncontrolled ion trajectories are not a significant concern because the 
background gas pressure is typically so high (le-3 < p < le-1 mbar typical) that the mean 
free path of ions in the IEC devices is less than twice the diameter of the device. Most 
ions in these systems are lost to the grid wires via scattering collisions or through charge 
exchange with the background gas [357]. The spatial representation of ions in these high 
pressure systems wind up resembling ion clouds, as opposed to well defined ion beams. 



Figure 5: OOPIC simulation of typical high pressure IEC discharge (p = 3e-3 mbar, Argon) 

Figure 5 shows a high-pressure (3e-3 mbar) discharge that is typical of most of the 
experimentation that is done [361]. The effect of scattering collisions is clear when 
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compared to the hard vacuum of figure 2. In this high pressure regime, the majority of 
fusion reactions are between ions and the neutral background gas. While the measurable 
neutron output is highest in this pressure regime, the scattering and charge exchange 
interactions of the ions with the background gas limits the ion confinement time. The 
maximum measured Q of IEC devices in this high pressure regime approached 10' 5 with 
deuterium and tritium (DT reaction) gas [18]. IEC operation in most experiments uses 
only deuterium (DD reaction). The average Q is of the order 10' 9 . 


1.2.2. OOPIC Modeling 

A Particle-In-Cell (PIC) modeling tool has been developed using a commercial-off-the- 
shelf (COTS) code called OOPIC. The OOPIC physics kernel was originally developed 
by UC Berkeley’s computational physics group [359, 360]. It is a 2-D code that has a 
wide range of built in functionality including a Monte Carlo Collision (MCC) model, 
simulated ion sources, realistic electron and ion bombardment ionization models for 
Argon, Helium and a number of other gases, and simple user definable geometries. 

PIC modeling informed much of the original investigations into improving ion 
confinement. It has been particularly useful for predicting ion bounce times, simulating 
the effects of multiple grids, and evaluating the development of the ion-ion counter- 
streaming instability in these low-pressure, non-neutral IEC systems with good particle 
confinement. 
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1.3. Scaling of UHV IEC Devices 


One of the particularly interesting characteristics of a non-neutral IEC system is the 
scaling. Contrary to the physical intuition of most power production systems where the 
power produced is proportional to the size of the device (for similar device types), a non- 
neutral IEC reactor has inverse scaling because achievable particle densities are limited 
by space charge effects. 


1.3.1. Derivation 

The physics of a non-neutral IEC device are dominated by the Poisson equation. 


V 2 ® = — Equation 1.1 

S Q 

We will define a dimensionless number from the Poisson equation that relates the 
maximum curvature due to particle space charge effects to the curvature of the vacuum 
potential structure which is just the Laplacian of the potential. We characterize the 
curvature of this Laplacian imposed by the boundary conditions on the grid wires with a 
typical radius, R typ and a typical potential difference, ® typ . We can then define our non- 
dimensional “IEC” number as follows: 


IEC = 


PchKp 

Sj&typ 


Equation 1.2 


It is convenient to make a number of assumptions for the purposes of comparison. We 
will first assume that R typ is the radius of the outermost grid (or the largest radius of the 
device), and we will assume ® typ to be the maximum potential difference between the 
anode and cathode grids in the device. We will also assume that any electrons generated 
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inside the IEC device will be lost on a time scale much faster than the ions, so the number 
density of electrons in the system is effectively 0. The charge density term can then be 
replaced by the maximum number density of ions in the system times Z e ff times the 
elementary charge, e. We will differentiate this number from the above number by using 
lower-case letters. 

eZ eff n i R] inode 

iec = — t 5525 r Equation 1.3 

8 o V® anode ~ ® cathode ) 

In order to maximize power density for a given system, it is always desirable to maximize 
the iec number, but space charge effects will limit the maximum achievable iec number 
in any non-neutral device. If we therefore assume that the iec number is constant for any 
system operating near the maximum space charge limit, the fundamental scaling 
relationships emerge. For two IEC devices with the same iec number, Z e ff, and grid 
potentials: 



Equation 1.4 


PowerDensity 



Equation 1.5 


Power 


Equation 1.6 


1.3.2. Numerical Validation of Scaling 

This very simple scaling argument stands up to test by detailed numerical simulation. 

The maximum ion density achievable in a system the size of our experiment (0.5 m 
diameter) at fusion-relevant cathode potentials (-100kV) is approximately 2e-14 m' 3 or an 
iec number of approximately 9. The simple scaling relationship indicates that in order to 
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achieve ion densities that are comparable to those in a Tokamak, the radius of the device 


must be decreased by a factor of 1000. 



Figure 6: Particle-In-Cell simulation showing ni 2.5e20m-3 at 0.5mm diameter device size 
When a geometrically similar OOPIC model was created with a diameter of 0.4 mm and 
the same -lOOkV potential was applied to the cathode, the peak ion densities were seen to 
top 2.5e20 m' 3 (an iec number of 7.4) as predicted by the simple scaling derivation in the 
previous section. 

Although the maximum iec number for a given size system does depend on the details of 
the design and implementation, for a system with good ion confinement that approaches 
its own space charge limits, the iec numbers are not expected to vary by more than one 
order of magnitude because fundamental space charge limits cannot be designed around 
in non-neutral systems. 
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2. Experiment Development 

Experimental validation of the hypothesis that ion confinement time in IEC devices could 
be improved by the addition of wire focusing grids was sought. It was decided at an early 
stage that this experimental validation would come from the simple geometric 
reconfiguration of the wire grids in a single experiment at high vacuum. The experiment 
would first be built in a conventional 2-grid geometry with a single cathode grid near the 
center of a much larger anode grid such as the geometries represented in figures 3 and 5. 
The ion confinement times for this system would be inferred from measured data. 
Additional focusing grids would then be added to the system (similar to figure 4), and ion 
confinement times would again be inferred from the new, “Multi-grid” data. These two 
data sets would then be compared to reveal whether they support the core hypothesis. 

A number of design and modeling techniques were employed to facilitate the design of 
the experiment. These models are detailed in the following chapter. The net result of the 
modeling effort was to confirm assertions that were inferred from the basic scaling given 
in chapter one and, more importantly, to quantify and predict the expected ion densities 
and confinement times of the system. 

The geometry of the experiment was largely dictated by experimental convenience - 
available equipment, measurement access, and ease of grid construction with available 
tools and techniques. This section discusses four critical aspects of the development of 
the experiment: confinement time detection techniques, sources of error, selection of 
primary diagnostic, and hardware development. 
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2. 1. Confinement Time Detection Techniques 


There are many ways to measure the confinement time of ions in an IEC system. This 
section will give an overview of the concepts that were considered for this experiment. 
All of the techniques presented require that the source of ions have either a well 
calibrated injection rate or a rapid on/off switching time. If neither of these conditions 
can be met, accurate measurement of the ion confinement time is not possible. 

2.1 .1 . Measurement of current to the grid wires 

As ions are injected into an IEC system where wire grids are used to establish the 
background potential structure, they will occasionally adopt a trajectory that will 
intercept one of the grid wires. When an ion impacts a grid wire it is rapidly decelerated 
and neutralized by the “sea” of electrons in the wire. In the process, the energy of the ion 
goes into heating the grid wire and potentially liberating secondary electrons from the 
surface of the wire. Because the ion is neutralized on the surface, the impact of an ion on 
a grid wire will draw an electric current which can be detected. The current drawn will 
be proportional to, but not identical to the ion current (defined by the rate of ion loss 
times the charge per ion). This discrepancy is due to the emission of secondary electrons 
from the surface of the wire which will cause the detected current to exceed the ion loss 
current by the average secondary electron yield factor. 


I grid = lions 0 + Tie ) Equation 2.1 

Early particle-in-cell (PIC) experiments showed that the number of ions lost to grid wires 
per unit time is roughly proportional to the number of ions trapped in the system at that 
time, particularly under high vacuum conditions where confinement limits are imposed 
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by collisions with background gas. Once the source of ions in turned off, this 
relationship will cause the rate of ion loss to follow an exponential decay curve. One 
method of estimating the confinement time of ions in the system is to numerically 
evaluate the time-constant of the assumed exponential decay. This can be done by 
measuring the time between the maximum loss rate and 1/e times that loss rate after the 
ion source has been turned off. In this manner, ion confinement time can be inferred 
from the measurement of the current to the grid wires as a function of time. 

2.1.2. Destructive Beam Dump 

Ion confinement time can be directly measured by a destructive count of the number of 
ions remaining in the trap a certain period of time after ion injection is turned off. This 
technique is referred to as a destructive beam dump (DBD). The DBD hardware would 
consist of an ion collector plate mounted in the beam path near the anode opposite the ion 
source. The collector plate would be held at or near ground so that ions could not impact 
the plate without substantial upscatter or thermalization of the distribution function. 

Then, when a measurement of the trapped ion population is desired, the potential of the 
plate would be rapidly lowered so that all of the ions would be energetically capable of 
impacting the collector plate. The location of the plate inline with the beam will cause all 
of the ions in the trap to be rapidly evacuated and neutralized on the surface of the 
collector. If the gradient of the electrostatic potential at the surface of the plate is such 
that it suppresses any net secondary electron emission, the time integral of the current 
that is drawn to the plate will be equal to the net charge that had been contained in the 
trap. When this measurement is made in conjunction with assumptions about the relative 
fraction of singly and doubly ions, the total population of ions in the trap at the dump 
time could be inferred with a high level of confidence. 
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DBD is appealing because it has the potential to be a direct measurement of the total 
trapped ion population at a given time. By taking many of these measurements with 
different delays from the termination of ion injection, a decay curve can be extracted of 
the actual ion population in the trap. Although many more measurements are required to 
perform this analysis, it is appealing in that the actual number of ions is measured as 
opposed to a signal that is simply proportional to the actual number of ions (as in the 
previous technique). 


The primary disadvantage of the DBD technique is that it is more expensive and time 
intensive to implement than the previous technique. In order to get a measurement of the 
ion population in the trap at a given time, the collector potential is rapidly lowered on the 
order of lkV for this experiment. Ideally, this change in potential would be 
instantaneous, but since the collector and its feed wires have a finite capacitance, an 
instantaneous change is not possible. The question then becomes, “how fast does it have 
to be?” This question can be answered by looking at the dynamics of the ions in the 
system and assuming something about the loss rate of ions from the trap. It is desirable 
to have the measurement of ions in the trap be “effectively instantaneous” - meaning that 
there is a negligible change in the trapped ion population during the time it takes to drop 
the collector potential from ground to the collecting potential. Mathematically speaking 
the condition is 


, N Trap 
' N Trnn 


Equation 2.2 
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PIC modeling of the system suggests that the confinement time of a 2-grid system is 
0(10) passes. Since it is desirable to compare the 2-grid system with the multi-grid 
system using the same diagnostic, the characteristic time for the collector potential drop 
( T collector drop ) should be less than half of the bounce time in order for accurate 
measurements of the 2 grid confinement time with a DBD technique. 



In addition, there is concern that the 2-grid beamline is not robust to minor asymmetries 
in the potential structure which could cause the beam to intersect the collector plate near 
the edge or miss the plate entirely. Clearly, if the plate is missed entirely, an accurate 
measurement cannot be made with this technique, but even if the beam impacts the 
collector near the edge, secondary electrons could be emitted from the collector to the 
ground which would introduce error into the “direct” measurement of trapped ions. 

There are other practical difficulties with implementing this type of detection scheme. 

The current required to change the potential of the plate due to its finite capacitance 
would need to be subtracted from the signal, and the need to rapidly pulse the potential of 
the collector over the very short dump time requires special power supplies. 

2.1.3. Capacitive Detection of Trapped Charge 

A probe or plate like the one used as a collector in the DBD technique can be similarly 
used to detect the trapped ions capacitively. Not only is this a direct measurement of the 
ions in the trap, it is also a non-destructive measurement so the entire decay can be 
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monitored with each injection without the need to take multiple measurements with 
different destructive delays on different injections as in the DBD. 

Capacitive detection of ions in the trap relies on the detection of the very small “image 
charge” currents that are induced in the probe or detector plate by nearby ions. In effect, 
the detector becomes one half of a capacitor and the trapped ions become the other half. 
The minor perturbation of the electric field in vacuum caused by the presence of nearby 
ions induces a small current in the wire leading to the plate if that wire is connected to 
ground (or some constant potential). Since there can be no electric field inside the 
conductor, the charges will rearrange to shield the field perturbation at the surface. That 
small current can be detected with an ammeter circuit. 

A primary difficulty in the use of a capacitive detector is that the measurement relies on 
an assumed location of the ion beam in order to estimate the effective capacitance of the 
detector. Based on the PIC models of the 2-grid device, the geometry of the beamline in 
this configuration is not at all well understood, so the effective “capacitance” of this 
detector is uncertain. It is therefore difficult to determine absolute magnitude of trapped 
charge using this technique, but temporal resolution should be excellent when using a 
small probe with <lpF capacitance. 

2.1.4. Laser Induced Fluorescence 

Laser induced fluorescence (LIF) of the trapped ion population could, in theory be used 
to measure the density of ions trapped in the system. LIF is appealing because it is a non- 
destructive direct measurement of the trapped ion density and it has the potential to 
directly show the spatial variation of the beam density. In LIF a laser beam is projected 
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through a measurement region that is anticipated to contain an ion population. The laser 
is tuned to excite a particular electron transition. The decay of that excited ion releases a 
photon in any direction. An aperture/filter arrangement in conjunction with a 
photomultiplier tube is then used to detect some small fraction of these photons. That 
signal is then assumed to be proportional to the density of ions in the region of space that 
is intersected by the beam and the field of view of the detector. By monitoring the density 
as a function of time (which may need to be done statistically with multiple 
measurements because of the small signal) it is theoretically possible to back out the 
trapped ion decay curve and hence the confinement time. 

Due to the expense associated with the development of a LIF detector, LIF was never 
seriously considered as a practical detection option for this experiment. It is mentioned 
here because it has substantial promise to reveal the shape of the actual beam envelope 
which could powerfully inform future investigations by validating or invalidating much 
of the PIC modeling that has been done. 

2.2. Sources of Error 

Every experiment has many potential sources of error which at best introduce inaccuracy 
or noise into the data and at worst cause the experiment to completely fail to measure 
what it is attempting to measure. This section details a number of the potential sources of 
error in what is anticipated to be levels of increasing significance. 

2.2.1. Earth’s Magnetic Field 

Charged particles will experience accelerations due to electric and magnetic fields. 
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Equation 2.4 


a = -^fe + vx2?l 
m 

While typical electric fields in the experiment are O(10 4 V/m), the v-cross-B term is only 
of order 10. This thousand- fold difference in the terms of the equation of motion is the 
justification for ignoring the effects of the Earth’s magnetic field in the PIC modeling. 
However, the effect is real and will result in the very slight broadening of an ion 
trajectory that would otherwise turn back directly upon itself in a linear trajectory. This 
skew in the trajectory that is introduced can be estimated by integrating the acceleration 
perpendicular to the nominal, radial ion trajectory over the period of one pass (or one half 
of the bounce time). This method implicitly assumes small angle deviations. 


Although in the actual experiment, the ion beam is out of alignment with Earth’s field 
lines by only 40 degrees (+/-10 degrees as measured with a compass on the device), in 
order to be conservative I will assume that the beamline is perpendicular to the local 
field. The average ion velocity is given approximately by: 



Equation 2.5 


and the total trajectory skew per pass is approximately given by: 
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Equation 2.6 


With a bounce time of 24 microseconds and a local field strength of approximately 55 
microTesla, this expression evaluates to a total skew of roughly 733 microns per pass for 
Ar 1+. Due to their faster average velocity, singly charged Helium ions, could be skewed 
by up to 2.3mm per pass, but these distances are arguably negligible since they are small 
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with respect to the grid openings at the anode (~18cm). It would require approximately 
123 passes (or about 2.9 milliseconds) for a singly charged argon ion to move from the 
center of the beam to an anode grid wire based solely upon this magnetic drift. A singly 
charged Helium ion would require 39 passes or approximately 147 microseconds to drift 
the same distance. Since confinement times have generally been limited to 
approximately 10 passes, the effect of this magnetic drift will be ignored. Asymmetries 
in the electric field are far more likely to cause significant perturbations of the ion 
trajectories than magnetic drift. 

While the impact on ion trajectories is negligible, the Earth’s magnetic field could 
substantially alter the trajectories of secondary electrons. But since all magnetic lines 
intersect the anode, even very low energy electrons with Larmor radii much smaller than 
the device scale will wind up impacting the anode. Due to the broad energy spectrum of 
secondary electrons, and the fast dynamics compared to the ions, electron confinement 
was not modeled. It is expected to be poor. 

In addition to directly impacting the particle trajectories in the experiment, the Earth’s 
field could also produce a Hall-effect current in the voltage divider resistors which could 
result in a local corona and flashover effect on the surface of the resistors. In order to 
protect the voltage dividers from corona and flashover and to prevent unwanted leak 
currents, the voltage dividers are potted in Silicon RTV. There is no evidence from the 
experiment that suggests that this phenomenon is occurring. 
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Overall, the Earth’s magnetic field is not expected to introduce substantial error into the 
measurement of the ion confinement time. 

2.2.2. Secondary Electron Emission 

In the commercial PIC code OOPIC, the Vaughn model is used to predict secondary 
electron emission from the surface of the grid wires [353]. There is large uncertainty, 
however, in how much of the energy of impacting ions results in secondary electrons 
versus surface heating versus sputtering. A thorough literature search was conducted, and 
the work of Szapiro and Rocca [355] was determined to be a good basis for our estimates 
of secondary electron production. 

Due to the possibility of a current cascade where secondary electrons from one grid wire 
are accelerated and impact other focusing wires at high energy which then produce their 
own secondary electrons, only the measurement of the current to the cathode grid will be 
used as being truly proportional to the ion flux. 

Since the cathode will be held at -lOkV for both the 2-grid and multi-grid cases, Szapiro 
and Rocca would estimate our secondary electron yield at roughly 3 electrons/ion for 
lOkeV argon ions impacting a stainless steel wire. Our stainless steel longitude wires are 
annealed in the atmosphere, however, there is an oxide layer covering the stainless steel 
wire which is likely to substantially enhance the secondary electron current. The 304 
stainless steel from which our wires are made is composed of Iron (66-74%), Chrome 
(18-20%), and Nickel (8-10%) with trace amounts of other elements. Based on the 
research of Allen, Dyke, Harris, and Morris [356] it is likely that the majority of the 
oxide on the wire surface is iron oxide. Direct data for argon ions impacting oxidized 
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iron at relevant energies were not located, but secondary electron yield from other oxides 
with lOkeV argon ions could be as high as 6 to 9 electrons/ion (for Aluminum oxide and 
magnesium oxide respectively) [355]. 

These yields are in stark contrast to the yields of secondary electrons off of clean, pure 
metal surfaces which are typically less than 1 electron/ion and depend substantially on 
ion incidence angle as well as energy and crystallographic orientation of the base metal 
[354], 

While there is a large amount of uncertainty in the actual yield of secondary electrons, the 
secondaries will actually increase the measured current to the grid in proportion to the ion 
flux (for the cathode). Therefore secondary emission will be simultaneously a source of 
increased signal strength and increased signal noise. Due to the statistical nature of the 
process and the high frequency of grid impacts and secondary emission, the signal 
strength is expected to increase much more than the noise. Of course, the oxide layer on 
the grids may sputter off which could cause the signal/noise to degrade with time as the 
grids are sputtered clean. The time of this variation will undoubtedly be a function of the 
depth of the oxide layer on the grids and the conditions at which the system is run 
(background pressure, applied voltage, injected beam current, etc.). 

2.2.3. Modeling Uncertainties 

Ions can be lost via undetectable pathways: recombination, or ion neutralization on the 
ion source itself. These loss mechanisms have not been modeled because they are 
assumed to be negligible compared to the primary loss mechanisms already discussed. If 
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confinement is worse than expected, it may be worthwhile to model these other possible 


loss mechanisms. 


2.2.4. Equipment Limitations 

The experiment described in this thesis was conducted in a vacuum chamber that was 
consistently capable of a high vacuum base pressure on the order of 5e-7 mbar. The 
lowest pressure achieved in the chamber after extensive cleaning and a 3 day pump-down 
was 8e-8 mbar. The facility does not have provisions for bake-out or an ion pump for 
UHV operation. While these base pressures should allow the demonstration of improved 
confinement in the multi-grid configuration, the chamber is not currently capable of 
achieving the base pressures identified [357] as necessary for ion-ion collision limited 
confinement. 


Grid Model 



Figure 7: Schematic model of one high-voltage grid circuit 


The four high voltage power supplies that control the grid potentials are also each limited 
to a maximum current of approximately 425 pA and a potential of -10 kV. The 
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maximum power output of each of the four independent supplies is approximately 4W. 
The geometry of the HV feed-through and the wire grids resulted in a built-in capacitance 
of approximately 3pF as shown in figure 7. Because 1/e ion confinement times could be 
as low as half the bounce time (or a minimum of around 3ps for Helium ions and a -lOkV 
cathode), the resistance between the grids and ground (R2+R3 in figure 7) could be at 
most 1MQ in order for the RC decay time constant to be lower than the decay time we 
are trying to measure using the direct measurement of current to the grid wires. This 
resistance limitation combined with the power/current limitations on the HV supplies 
resulted in a minimum cathode voltage of only -400V for the grid-current diagnostic. 

It should be noted that by floating (isolating) the potential of either the power supply 
itself or a difference amplifier across R1 in figure 7 the current could be measured 
without this restriction on the grid voltage. Due to practical constraints on isolating these 
components in the existing laboratory environment, these options were not pursued. 

The use of the destructive beam-dump diagnostic was prevented by the unavailability of a 
second functional fast-switching HV supply. One is necessary to control the injection of 
ions, and another would be required to dump the trap after a certain amount of time. 

Only one functional fast-switch was available for the present investigations. A second 
switch was on hand, but it stopped working before useful measurements could be made. 

Differential pumping was not available, so it was not practical to maintain a substantial 
pressure differential between the ion production region and the bulk of the chamber 
volume. Therefore, the strength (current out) of the ion source is proportional to the 
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background gas pressure in the chamber. When the pressure is reduced to levels where 
multiple order-of-magnitude improvements in 1/e confinement time are predicted by the 
particle in cell models, the ion source is so weak that the signal to noise ratio is too low to 
discern the 1/e confinement time from the data. 

A final, and potentially quite important source of error arises from the construction of the 
ion source itself. The ion source is contains a thoriated tungsten filament which serves as 
an electron source. The filament is located outside the anode grid. A fine (~lcm 
spacing) wire mesh screen covers the equatorial section of the anode. See figure 18 in 
the next section for a schematic drawing of the ion source. The filament potential is 
rapidly (<1 ps) switched from the ground (anode) potential to -150V. Electrons emitted 
by the filament are then accelerated toward the anode screen. Because the screen is 87% 
transparent, the Hirsch formula [21] suggests that electrons will make 3.58 passes back 
and forth through the screen before impacting one of the wires. Background gas in the 
region is ionized via electron bombardment. Ions that are generated inside the anode grid 
fall into the potential well of the device, while ions that are created outside the anode grid 
are accelerated towards the filament and are expected to generally result in the production 
of more secondary electrons. When the ion source is turned off the filament potential is 
rapidly (<1 ps) raised to ground. The electrons that are now “boiled off’ the surface 
thermionically do not have sufficient energy to ionize the background gas. The ion 
source is effectively turned off. The possible source of error arises from ions that were 
generated on the outside of the anode just before the filament potential was grounded. 
Those ions may be extracted from their location outside the anode screen into the 
potential well by field penetration from the inner grids through the screen. This situation 
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could result in a “slug” of ions born outside the anode entering the well after the filament 
potential is raised to ground. This transient slug of ions could result in an increase in the 
measured 1/e confinement time. There is evidence to suggest that this happens in the 
experiment. 


2.3. Selection of Primary Diagnostic 

Although it was originally anticipated that the measurement of the current to the grid 
wires would be the primary diagnostic for the experiment, the limitations described in the 
previous section resulted in a change in primary diagnostic. 


2.3.1. Diagnostic Comparison 

The candidates for the primary diagnostic were: the measurement of the current to the 
grid wires, the measurement of the current to a probe near the anode held at a potential 
lower than the anode, and the capacitive pickup on the probe near the anode held at a 
potential above the anode. As previously mentioned, DBD and LIF had been discounted 
due to equipment limitations. 

The measurement of the current to the grid wires is limited to cathode voltages of -400 V. 
There was concern that with this relatively high cathode voltage and correspondingly 
shallow potential well, electrons born on the filament at -150V could penetrate quite far 
down into the potential well. The ions born in the well from the electron bombardment 
of the background gas would then have a relatively broad energy spectrum compared to 
the well depth. This sort of “spread ion source” was shown in OOPIC simulation to 
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result in less well confined ions. The birth location and potential of the ions was 
observed to have an impact on the ion confinement time in the OOPIC simulation, so this 
detection technique was not selected as a primary diagnostic. It is useful for purposes of 
comparison however. 

A probe inserted into the ion beam opposite the ion injector was used to make two other 
types of diagnostics (see figure 19). When the probe is biased at a potential lower than - 
150V, all ions in the system are energetically capable of being neutralized on the surface 
of the probe. The probe at this potential can act like a separate “collector grid” and allow 
the direct measurement of lost ion current similar to the measurement of the current to the 
grid wires. This type of diagnostic is appealing due to the ability to disconnect the grid 
voltages from the collector voltage. The cathode grid can be biased to a potential much 
lower than the -400V limit of the previous diagnostic, and the time resolution on the 
detector can be improved by the use of a diagnostic with a smaller inherent capacitance 
(<lpF) and a lower resistance (<1MQ). The very presence of the probe, however, 
significantly perturbs the potential structure around it because the Debye length of the 
system at the operating densities is significantly larger than the diameter of the device. 
This situation is true for any non-neutral IEC device. The perturbation of the potential 
field around the probe results in an asymmetry which has the ability to change the 
confinement properties. The biggest disadvantage of the use of the probe in this manner 
is that it introduces an intentional ion sink into the ion beam that would not exist in an 
operating device. Due to the 2-D nature of the OOPIC code, the effect of the probe at 
this potential could not be realistically modeled because in the simulation, every ion that 
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approached the probe would neutralize on the surface. While this scenario is not 


physical, there is not a good way of quantifying the error with the existing tools. 

When the same probe is biased slightly above the anode potential (+0.3V, figure 20), ions 
bom in the system are energetically incapable of impacting the probe. The ion sink is 
removed. A significant asymmetry in the potential structure still exists, but OOPIC 
modeling (which could be used in this case since no ions are lost to the probe) suggested 
that a potential structure with good confinement would still tend to have good 
confinement with the probe in the system on one side. The probe at this potential acts as 
a capacitive pickup. Although the currents to the probe will be lower, the probe should 
not significantly reduce the confinement time of ions in the system. This diagnostic was 
therefore chosen as the primary diagnostic for these experiments. 


2.4. Hardware Development 

Significant effort went into the development of the facilities necessary to show the 
potential for multiple grids to improve the ion confinement time in IEC systems. This 
section details some of the important work should it need to be recreated in the future. 

2.4.1. Vacuum Chamber Retrofitting 

A stainless steel vacuum tank was acquired from another university where it had been 
used for sputtering materials. The inside of the tank was blasted clean and a hinged door 
was retrofitted. An additional port was added to the tank to accommodate the 1000 L/s 
turbopump used for all high vacuum operation. A leaded-glass window was added to one 
of the ports. A lOkY electrical feedthrough was added to supply the grids, along with a 
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combined gas and electrical feedthrough structure for the ion source and a 700V 
diagnostic feedthrough opposite the ion source. Finally, a stand was constructed that 
would allow the chamber to be transported and supported horizontally for easy access. 

Figure 8 shows the chamber in its final configuration. A combined pirani and ion gage 
pressure sensor was mounted to another port on the chamber allowing pressure 
measurement from 1 atmosphere to 10' 10 mbar. 

Power supplies and diagnostic hardware are visible in the rack on the right. A 2 GHz 
Tektronix 2014 Digital Oscilloscope was used for data acquisition. Output was through 
an RS-232 serial cable to the computer (monitor visible on the bench at the far right). 



Figure 8: Vacuum chamber 
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2.4.2. Grid Fabrication 



Figure 9: Fabrication of grid wire jig 


Since the decision had been made to construct the multi-grid portion of the experiment 
with five concentric spherical wire grids, it was necessary to develop a standard 
methodology for constructing the grids. The final solution that was arrived at is reported 
in this section. 


After considering many different options for grid wire materials, 1mm diameter 304 
stainless steel wire was chosen for all grids so the structures would have sufficient 
rigidity to survive normal handling, and so conventional, spot-welding techniques could 
be implemented with ease. It is expected that a more detailed analysis of the material 
choice will be required for future devices. 
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A jig to hold the 304 stainless wire in place at the proper curvature for spot welding was 
machined from a piece of polycarbonate as shown in figure 9. The wire was laid in the 
trough and spot-welded to the desired diameters. It was then removed from the jig. The 
wires that were to become longitude lines were then annealed so they could be cut and 
still maintain the desired radius of curvature. After some limited experimentation shown 
in figure 10, it was discovered that running a current of 35 A through the hoops for 1 


mi nute would give the desired effect. 



When less current was used for less time, the wire would show some residual “spring- 
back,” but if more current was used for much more time, the wire would tend to deform 
significantly under its own weight and the forces applied by the alligator clips used to 
supply the current. 
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Figure 11: Atmospheric annealing of longitude wires 
Figure 1 1 shows the atmospheric annealing of one of the grid wires. Once the annealing 
of the longitude lines was complete, they could be cut in sections and maintain the 
desired radius of curvature. The latitude lines were not annealed because they did not 
need to be cut during the assembly procedure. The longitude lines had to be cut so the 
feed-throughs supplying the inner grids could pass through the pole and introduce the 
minimum possible field asymmetry. 

In order to attach the longitude wires to the polar feed-through, special stainless steel 
rings were fabricated from 1/16” plate using an OMAX waterjet cutter. These rings 
provided structural support and an electrical connection for the grids. 
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Figure 12: Spot welding longitude wires to retaining rings 
The eight longitude wires were spot-welded to the support rings as shown in figure 12. 
The support rings were designed to fit onto custom machined aluminum parts that were 
held onto the alumina feed-through stalk at the appropriate radius by a set screw. 

After the longitude lines were welded to the retainer rings, the latitude lines would be 
spot welded onto the outside of the longitude lines at the appropriate position. A simple 
wire spacer was used to locate the hoops concentrically around the poles. Errors in the 
actual location of the longitude lines could result in a slight perturbation of the latitude 
line position. In general, this construction error was measured to be less than 5%. 
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Figure 13: All finished grids 


Figure 13 shows all of the completed wire grids. Slight asymmetries in the grid wires are 
visible to the naked eye, but position errors of the latitude line locations were all 
measured to be less than 8% of the respective radius. Longitude line errors were 
significantly less than that. 

In order to install the grids in the final multi-grid system, they had to nest one inside the 
other. This was accomplished by again cutting the longitude lines near the equator this 
time. The cut was made actually just “south” of the equator - right above the southern 
“tropic” latitude line. This asymmetric cut was used to improve the ease of re-assembly. 
It was desired to have one well-located, stiff wire end, and one more easily manipulated, 
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less well-located wire end that would be moved to align with the stiff end. This was 


assumed to be easier than lining up to wire ends of moderate stiffness. 



Figure 14: The "Awesome Blossom" 

The resulting “half-grids” resemble a flower blossom as shown in figure 14 and a radar 
dish or directed listening device as shown in figure 15. 



Figure 15: The "radar" with alumina feed-through and aluminum ring supports 


45 


When the grids are assembled in the chamber, 1” sections of 1/16” diameter stainless 
steel tubing are crimped onto the longitude wire ends on the “radar” half, which allow the 
insertion and perfect alignment of the longitude lines on the upper, “blossom” half. 

Figure 16 shows the completed multi-grid IEC assembly inside the vacuum chamber. 



Figure 16: Multi-grid IEC in the vacuum chamber 


2.4.3. Ion Source 

During the course of these investigations a number of different ion source geometries 
were tried. The final ion source used for all of the experiments reported in this thesis 
was designed to have a minimal impact on the potential structure inside the anode grid. 


46 




Figure 17: Ion source used in experiment 

Figure 17 shows a picture of the ion source used for all experiments reported herein. Ions 
are generated by electron bombardment of the background gas in the chamber. The 
electron source is a thoriated tungsten wire filament located outside the anode grid. The 
filament continually emits electrons thermionically. The filament is powered by an 
isolated variac that is always on. When ion generation is desired, the filament potential is 
reduced from ground to -150 V by a fast amplifier capable of driving the filament circuit 
to the desired potential within 3 ps. Electrons emitted from the filament at the lower 
potential are accelerated towards an 87% transparent screen mesh which is held at 
ground. Figure 18 shows a schematic of this ion source. 
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The peak of the electron bombardment ionization cross-sections for both Argon and 


Helium lie in the range of 100-150eV electron energy. Ions are therefore generated all 
around the screen mesh which is positioned such that it is tangent to the surface of the 
spherical anode grid. Those ions which are generated inside the anode grid fall into the 
well, while those generated outside are accelerated towards the filament and aluminum 
foil “neutralizer” which is electrically connected to one of the filament leads. 



Figure 18: Schematic representation of ion source 

This design results in the quick removal of ions generated outside the anode while 
allowing the ions generated inside the anode to fall unobstructed into the IEC potential 
well. 
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One drawback of this source is that after the filament potential is returned to ground, 
there are ions that were generated between the filament and the screen which now are not 
attracted to the neutralizer and are slowly extracted into the IEC potential well by field 
penetration through the grid mesh. Due to the high electron density in this region, there 
exists a large number of these ions relative to the ions generated on the other side of the 
screen. This results in a “slug” of ions at shutdown which can be seen clearly in the data 
when the potential structure has good confinement. 

2.4.4. Ion Detector 

The capacitive probe was the primary ion detector used in these experiments, although 
measurements of the current directly to the grid wires were also made at grid potentials 
near ground (down to -400 Y). Figure 19 shows a picture of the capacitive probe wire 
poking through the center of the anode mesh. 
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The detector circuit is a simple 1/100 voltage divider on a 30V DC power supply which 


held the probe at a no mi nal potential of around 0.3V (actually 0.28 V due to resistor 
tolerances). A schematic of the detector is shown in figure 20. 

Since the maximum birth potential of any ion in the system is 0V, this 0.3V offset would 
ensure that a negligible number of ions will be neutralized on the probe due to the 
thermal spread of the ion energy at birth (assumed ~.03eV). Most ions should be bom at 
lower potentials (an average of around -25V based on the Argon cross-section). Of 
course, collective effects such as the two-stream instability have the potential to spread 
the ion distribution function on a much faster time-scale than the collisional 
thermalization time which could result in some ions impacting the probe, but due to its 
small area compared with the anode grid and screen, it is highly unlikely that the probe 
will be a significant ion sink. 


I EC potential well 



Ion turn-around 



Scope Thevanin equivalent 
C = 95 pF 
R = 1 MQ 


-±r~ 30VDC 


Figure 20: Schematic of capacitive detection circuit 
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The Tektronix TDS 2014 scope was attached to the wire probe on the middle of the 
divider via a IX scope probe. With these resistors, the estimated maximum bandwidth of 
this detector is around 100kHz which is about at the limit for Argon (with a bounce time 
of ~23 ps). Lower value resistors were used to search for the predicted two-stream 
instability in Helium due to Helium’s 8 ps bounce time, but lower value resistors reduce 
the strength of the signal as well as increase the bandwidth, and Helium had a much 
weaker signal to begin with due to its lower ionization cross-section. 

The following chapter presents some of the computational modeling that was done prior 
to and concurrent with the design and operation of the experiment. 
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3. Computational Modeling for Design 


In order to estimate IEC device performance, two independent computer codes were used 
to predict ion behavior. The commercial PIC code called OOPIC was introduced in 
chapter 1 . In addition, a custom code was written in MATLAB to approximate the true, 
3-D potential structure. The potential variation along the center of the beamline was then 
used in conjunction with the 1-D paraxial ray approximation to solve for the beam 
envelope and the maximum confined core ion density. This chapter explains these 
computational tools in greater depth. 


3. 1. Semi-analytic potential model for a multi-grid IEC device 

A means of solving for the full, three-dimensional potential structure around a multi-grid 
IEC device was desired. Conventional finite difference approaches are not convenient. 
The small wire size and the large volume and wire spacing make for a computationally 
inefficient solution. Finite element methods with adaptive grid spacing were not 
seriously considered due to the scope of the computational problem. Instead, a semi- 
analytic approach was adopted with a truncated series potential representation. 

3.1.1. Derivation 

First the space charge of the trapped ions in neglected, so the potential structure is 
assumed to be a solution to the Laplace equation. 


Equation 3.1 
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We choose to represent the potential in spherical coordinates where 0 = longitude and 
</> = co - latitude . The potential and its derivatives can then be represented by a 
summation of orthogonal spherical harmonics 

® T"! U' km r k + (cos^)cos m0 + F t Equation 3.2 

k=Om=-k 

and 

^^-(r,^,#) = _ (A+2) ]^T (cos cosot ^ Equation 3.3 

5/* k=Om=-k 

where P™ is the associated Legendre function of the first kind and the index is an 
integer from 1 to n+l where n = # of concentric spherical wire grids. The coefficients A 
and B (and the constant offset F) uniquely determine the potential. It can be inferred 
from equation 3.2 that for a finite potential solution, all A n+1 = 0, and similarly, all B ! = 0. 

A finite potential solution also requires A — > 0 and B — » 0 as k — > oo . Therefore a 
truncated series with sufficiently large ‘T’ should closely approximate the exact solution. 
The problem becomes one of solving for the A and B coefficients for the finite series 
given the appropriate boundary conditions. With this type of representation, the 
boundary conditions become trivial to impose. From continuity we get: 

F i+l = F i Equation 3.4 

B l k m ~ B’km = (a^J — A‘ km Equation 3.5 

where r, is the radius of the /th grid. We then define an angular delta function such that 

rJS(0-jS)d0 = 1, (0)8(0- P)d0 = l^-,0<p<a Equations 

0 0 r i 

Poisson’s equation (the flux condition) at the grid wires can then be expressed as 
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jn 

t + 1 


) + Y d q t S{e 



Equation 3.7 


where q t is the charge per unit length on the grid wires, t = # of latitude wires , and 
1 = # of longitude wires . Note the form of this equation assumes an even angular 
distribution of both longitude and latitude lines. From this point forward, the expression 
on the right hand side of equation 3.7 will be referred to simply as RHS and the 

50 

expression on the left will be referred to as LHS. Substituting — on the LHS we get 


Z Z [ kr i (k ' 1) ( A L - A L')-( k + 1 K (k+2) (B‘ m - B‘J )]/>* (cos </>) cos mO = RHS Equation 3.8 


Then using continuity as expressed in equation 3.5, this expression can be simplified to 


^ ^ (2k + 1 )r j (k l) \A' km - A k J )P k m (cos (f>) cos m6 = RHS Equation 3.9 


or 

y 5^, - (2k + 1 )r ; ' (k+2) [B' km - B‘ kn ] jP k ‘ (cos <f>) cos md = RHS Equation 3. 10 

It can be clearly seen that the RHS expression uniquely determines the relationship 
between the A 1 and A 1 ' 1 coefficients, and similarly B t+1 and B 1 coefficients. Because we 
know that all the B' km and A km 1 coefficients must be zero in order for the series to 
converge at r = 0 and r = oo respectively, all of the A and B coefficients can be solved 
based upon the RHS expression. These coefficients can be extracted from the 
summations because of their respective orthogonality. In particular, cosine orthogonality 
can be used to simplify the expression. This step will eliminate the sum over m. 
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| (LHS) cos m'0d0= \{RHS)ca&m' OdO 


Equation 3.11 


This picks out m=m ’ and gives a factor of pi. Similarly, the orthogonality of Legendre 
polynomials can be used to eliminate the sum over k. When complete, equation 3.9 
becomes 

27ir < i k ~ l) (k + m )- ^ 4 ^ _ jr- 1 ) = J | ( RHS)Pf (cos (f)) sin (f > cos m' 6d Od(f) Equation 3.12 

C k-m)\ 0 0 

Substituting in the RHS term and integrating reveals the simplified relationship between 
the coefficients. 

2 ^r/k-') + - A'J ) = — — S m0 2 kY P™ (cos yfr) + S m(zl) l f P” (cos </>) sin 

( k-m)\ £ 0 r i L * 

Equation 3.13 

f 1, a = b 

where 8 ab = \ and z is any positive integer. So to solve for any coefficient, 

[0, a * b 

simply start with A ^ = B[ m = 0 and use. 


km 2 ne 0 r* ( k + m)\ 


feH 2 7rs 0 r- (k+l) (k + m)\ 


where X = <5 m0 2 tt^ P™ (cos S m(zl) l j P™ (cos ^)sin (f 


Equation 3.14 


Equation 3.15 
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In this manner the vacuum potential structure around an arbitrary number of 
independently charged concentric, spherical wire grids can be modeled in three- 
dimensions. 


The charge q t on the wires is approximated from a simple spherical capacitive 
relationship multiplied by a factor to account for the high transparency of the “spheres” 
and the distribution of charge over the surface area (proportional to r 2 ). 


q, = 4 ns c 


%A (®, -<D, +1 ) 


1 _ 

r. r . , 


Equation 3.16 


In the coded implementation of this model, Gaussian distributions of charge in (/> and 6 
were used instead of delta functions to limit the bandwidth of the spatial representation of 
the wires to allow for faster convergence with a smaller truncation value of k. More 
information about the details of the modeling can be found in the code contained in the 
appendix. 


Because a representation of the potential in all three dimensions is difficult to plot, 
sample cross-sections are presented in this chapter either at constant r, (j), or 6. In general 
a kmax of 32 yielded smooth partial derivatives in r which were required for paraxial ray 
beam envelope approximations. Figures 21, 22 and 23 show the modeling results for a 
typical potential variation on a plane through the equator. 
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Figure 22: Different view of same potential map 


Figure 21: Potential map of equatorial cross-section (§=n!2) of a multi-grid IEC device. 
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Figure 23: Same potential model, "pie" section from beamline to grid wires 

3.1.2. Limitations and Inaccuracies 

As with any model, there are limitations to the range of applicability of the model 
described in the previous section. In particular, the model has two fundamental 
inaccuracies which impose strict limitations on its range of validity: 

1) Truncation of the infinite series results in unrealistic high frequency “noise” in 
the potential representation - especially when approximating field gradients 
over very short distances without applying any type of spatial potential 
smoothing near the cut-off frequency. 

2) The assumption of constant charge density on the grid wires in wrong. In fact, 
the charge density varies significantly from the equator to the poles and in 
regions where wires overlap. The proper boundary condition is a constant 
potential on the grid wires. Unfortunately, that type of B.C. cannot easily be 
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implemented with this modeling approach. Plots of equipotential surfaces in 
figure 24 allow assessment of the significance of this inaccuracy. 




Figure 24: Spatial maps of the potential at a constant radius in between latitude and longitude wires 
on the equator, note unrealistic high frequency noise in the equipotential representation and 
potential deviation at wire overlap at the corners of the figures. 
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Figure 24 illustrates these inaccuracies clearly. An exact model would show no 
difference in the potential between the corners and the middle of the wires (along the 
edges of the figures). An exact model would also have smooth equipotential surfaces 
instead of the ripples introduced by the truncated series approximation. 

These inaccuracies, together with the assumptions in the model itself result in the 
following limitations: 

1) If resolution of the potential and derivatives of the potential are required at a 
given spatial frequency, the chosen value of kmax must not only satisfy the 
Nyquist condition, it must be twice the typical Nyquist criterion (at least 4 
times the frequency) to allow the spatially smoothed signal to satisfy the 
Nyquist criterion i.e. the smoothing window should correspond to Vi the size 
of the minimum detail desired, and kmax must be Vi the wavelength of the 
smoothing window. 

2) The potential model should not be used at all within the minimum resolution 
(as defined in #1 above) of the grid wire locations. 

3) The potential model should not be used within the range of the inaccurate 
boundary conditions. This range of inaccuracy can be seen in the 
equipotential plots. When the equipotentials do not resemble the bulk shape 
of the grid wires (e.g. near the wire overlaps), the model will not give an 
accurate potential. 

Fortunately, this model is to be used only for the purposes of preliminary experiment 
design, and for that reason, the only potential that is of significance is that along the 
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center of the beam line. This model should yield reasonable estimates of the potential for 


these crude purposes. 




Figure 25: Constant 0 cross section plots of potential structure of a multi-grid IEC device on a grid 
wire v. on the beam path 

Figure 25 illustrates the continuous potential (upper plot) and derivative of the potential 
(lower plot) along the beam line and the discontinuous derivative of the potential at the 
grid wire locations as a result of the charge located there. The upper plot also shows the 
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regions where the potential structure is “confining” to ions and where it is “unconfining” 


i.e. near the lowest potential grid. 


3.2. 1-D Paraxial Ray Equation Approximation 

In order to estimate the maximum confined non-neutral ion density achievable in a multi- 
grid IEC device, a simple model of the variation of the radius of an ion beam based on 
the paraxial ray equation was developed. 


3.2.1. Model Description 

The output from the semi-analytic potential model described in section 2. 1 was used to 
approximate the potential and the derivatives of the potential along an equatorial beam 
path in the multi-grid IEC device. The paraxial ray equation was used to predict the 
evolution of the ion beam envelope based upon that variation in the axial potential. This 
model was derived primarily from the work of Humphries [362], Since there are no 
magnetic fields in an IEC device, the general form of the paraxial ray equation reduces to 


®'R' _ 0"R s 2 K 
20 40 R 3 R 


Equation 3.17 


where R is the radius of the beam envelope (perpendicular to the beam line and the r 
vector in the IEC device), sin this equation is the beam emittance, and K is the 
generalized perveance defined relativistically as 


K = 


el 0 

2 7r£ 0 m 0 (Pycf 


Equation 3.18 
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In this equation, y = . , c is the speed of light and (3=v/c the ratio of the ion 

velocity to the velocity of light. The generalized perveance is proportional to the number 
density of ions, and it is the term that accounts for the space charge of the beam in this 
simple model. Although it is inaccurate, zero emittance (the laminar beam 
approximation) is assumed because the goal of this model is to predict the maximum 
space-charge limited confined ion density in the core of the device, not to predict a 
particular beam envelope. 

3.2.2. Numerical Integration of Beam Envelope 

Symmetry of the ion beam in the core of the device (a zero slope condition) is assumed 
along with a core radius. These assumptions provide the initial boundary conditions from 
which the paraxial ray equation can be integrated. MATLAB’s “ode45” routine is then 
used to integrate the paraxial ray equation from the core out to the radius of the injector. 

In order to verify convergence of the integration, the numerical integration is repeated in 
reverse, this time the initial conditions are the final conditions of the 1 st integration and 
the beam envelope is integrated on the way in. The deviation of the final radius and slope 
of the beam envelope is compared to the assumed initial radius to verify convergence of 
the routine. 

Figure 26 shows an example of a beam envelop that did not converge, while figure 27 
shows a converged beam envelope. 
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- Outward Beam Envelope for n c =1.37E-K)14 m-3, r Q =1.1 IE-002 m 

- Inward Beam Envelope for n Q =6.04E-K)15 m-3 t r Q =3.61E-002 m 



Figure 26: Improperly converged integration example 

3.2.3. Iterative Solution to Maximum Beam Density 

A maximum acceptance window which defines how close the modeled beam is allowed 
to come to the grid wires is assumed (typically 70-90% of the height of the grid wire 
from the beam axis). The initial beam current (ion density) is assumed to be two orders 
of magnitude lower than the density at which space charge effects are expected to be seen 
based on PIC modeling and simple back-of-the-envelope Child-Langmuir calculations for 
this size device. The beam envelop is then solved out to the injector radius. Near the 
injector radius the paraxial ray equation loses all validity because the implicit 
linearization assumes beam-like behavior which is not the case near the ion injection 
potential. If the beam envelope does not exceed the acceptance window at the grid-wire 
locations, the core density is increased by a factor of two, and the integration is repeated. 
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This iterative process is repeated until the space charge effects of the beam force the 
beam envelope outside the acceptance window of one of the grids at which point the 
routine is stopped. The maximum confinable ion density is thereby found to within a 
factor of two. 




Figure 27: Potential map (top) and beam envelope (bottom) versus beam-line position for very small 
core focus and high core ion density (not reproducible in OOPIC) 

Given the assumptions implicit in this method, it was expected to yield results only 
accurate to within an order of magnitude of the actual space charge limit, and when 
exceedingly small core radii were assumed, the maximum predicted density differed from 
the PIC model by more than an order of magnitude. When the assumed core radius was a 
significant percentage of the innermost grid radius, comparison to OOPIC modeling 
revealed agreement to within a factor of 2. Figures 28 and 29 show the similarities in 
the predicted beam envelopes from the modeling described above and the OOPIC 
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simulation. It is worth noting that OOPIC has only a 2-dimensional model of the 
potential, so minor deviations between the two codes should be expected. 


♦ Injection Point 

Beam Envelope (ot n c =7.50E+O13 m-3, r o =!.00E-002 m 



Figure 28: Model with large assumed core size that matches OOPIC results to within a factor of 2 
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3.3. Simulated Annealing IEC Design 

Due to the large number of independent variables involved in the design of a multi-grid 
IEC experiment, a code was developed based upon the beam modeling in the prior 
section to help find an optimal or near optimal design for a multi-grid IEC reactor given 
the limited resources of our group. Simulated Annealing (SA) is a design methodology 
for complex systems which depend on multiple independent design variables. SA is an 
active area of research in complex systems, and it has been shown to be an effective 
optimization tool for the design of complicated systems. Because of the large number 
of independent design variables (grid positions and potentials) in a multi-grid IEC, 
simulated annealing design was explored. 

The premise of simulated annealing design is that while there can be a very large number 
of independent design variables, the quality of any design can be quantified by a singular 
metric known as an energy function. The lower the energy function, the better the 
design. Because the overall design space can be very complex, there can be many local 
minima which could fool a gradient-search based design code. An SA approach initially 
resembles a Monte-Carlo method in that design variables are changed at random and 
even if a given change results in a higher value of the energy function, at the beginning 
that change is kept and another change to the design is made. The best design found is 
always kept in memory. As the code progresses, the bounds of whether a new design 
change is kept are gradually reduced until the code simulates a traditional gradient search 
method. This approach which gradually transforms itself from a Monte-Carlo design tool 
to a gradient search method has been shown to be quite good at finding designs at or near 
the global minimum energy function. The gradual decrease in the accepted bounds of 
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energy resembles the physical process of cooling (hence “simulated annealing”) - if it is 
done sufficiently slowly, the global minimum energy (single crystal) state will emerge. 

3.3.1. Energy Function 

For the multi-grid IEC design problem, an energy function was desired which would give 
a design with the maximum fusion energy output given the constraints on the design 
variables. The energy function was therefore set as proportional to the inverse of the 
estimated fusion power output from the device. 

EnergyFunction = — — 7 ^ — - — Equation 3.19 

n co re \®inj -^corefKore 

The fusion reaction rate is proportional to the square of the ion density in the core and the 
cube of the core size. The key implicit assumption here is that the fusion cross section is 
proportional to the square of the core ion energy, this assumption is only valid up to a 
certain fraction of the energy of the maximum fusion cross section for a given reaction 
(~50keV for deuterium). Since the grids in our experiment are limited to -lOkV 
potentials by the power supplies available, this assumption should be good for the current 
application. If future experiments are to use the same code, an energy function which 
incorporates the actual variation of the fusion cross-section for the desired fuel would be 
better. 


3.3.2. Design Variables and Limitations 

The radius and potential of each grid are all independent design variables. In addition, 
the radial location of the injector could be moved. The anode radius and the number of 
latitude and longitude wires were fixed, and the core radius was assumed to be a constant 
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fraction of the innermost grid radius. The innermost grid radius was limited to 10 cm 
maximum, and the grids were not allowed to be positioned within 2 cm radially of one 
another. These constraints were imposed for ease of fabrication, diagnostic access, and 
assembly. 


3.3.3. Summary of SA Results 

The simulated annealing design code converged on an optimal solution after 1600 
independent and random design perturbations. 


SA convergence history 



Figure 30: Final Convergence History for Simulated Annealing Design 
The initial design guess had an energy function that was nearly an order of magnitude 
higher than the final solution. Figure 30 shows the progression of the design with each 
successive iterative perturbation. The code successfully resulted in a grid design and a 
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recommended potential structure that would result in a significantly higher fusion power 


output than the initial guess. 



Figure 31: Number of occurrences v. energy in final SA run 
The progression of the “cooling” of the design can be viewed either by plotting the 
energy v. the iteration number, or by looking at the number of occurrences at different 
energy levels and the average energy level for all iterations at a particular pseudo- 
temperature as shown in figure 3 1 . 


The “freezing” of the design can be seen by the drop in the “entropy” and the rapid rise in 
the “specific heat” of the design as illustrated in figure 32. 
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Simulated Annealing Evolution 




Figure 32: Evolution of the SA design showing the freezing of the configuration 
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The final output of the S A code gave the potential structure and beam envelope shown in 


figure 33. 



Modified SA Results 



Figure 33: Final Output of SA design showing largest allowed inner grid and bunching of 
accelerating grids at coded limit of 2 cm. 

As can be seen in figure 33 above, the SA design code pushed the inner-most grid out to 
the coded limit of 10 cm. It also bunched the inner grid and the two accelerating grids to 
the minimum coded spacing limit of 2 cm. These can be explained by simply 
maximizing the size of the fusing volume and minimizing the spacing between the 
accelerating grids to reduce the space charge effects by maximizing the local field 
curvature near the accelerating grids. The injector potential is very close to ground, and 
the core potential is very close to the minimum of -lOkV. In fact the second and third 
grids out from the center are at - lOkV and 0V respectively - the maximum possible 
potential difference and the minimum spacing. It is interesting to note that neither the 
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first grid nor the fourth grid have any significant role to play in shaping the potential 


structure. This suggests that the optimum number of grids may be 3 (not 5). It is 
interesting to note that the while the energy function did change substantially, little of 
that change was due to a difference in density. In fact, the energy function showed the 
most sensitivity to the size of the innermost grid (and hence the assumed core size) and 
the difference between the core potential and the injection potential. The maximum 
density for Argon ions in the system was always seen to be around 8el3 m' 3 . 

While this tool may be useful for future investigations of actual fusion reactors, it was 
less useful than anticipated for informing the design of the present experiment. This was 
partly due to relating the energy function to possible fusion output instead of relating it to 
a measurement of the confinement experiment. Because most of the signals in the 
experiment will be proportional to the number of trapped ions, and that number does not 
seem to vary substantially given the maximum potential difference and fixed dimensions 
desired for experimental access, the actual hardware design was based off of 
experimental convenience instead of the results of this analysis. The code is included in 
the appendix of this thesis for future reference. 


3.4. Modeling with commercial OOPIC Pro code 

Throughout these investigations, extensive use of the OOPIC Pro code was made in order 
to gain insight into the detailed physics of the multi-grid IEC device. OOPIC Pro is a 2D 
planar particle in cell model with Monte Carlo collisions [358]. This section details some 
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of the predictions that were made prior to experiment operation based upon this useful 
tool. 

3.4.1. 2-Grid Modeling and Confinement Estimations 

Initial computational investigations explored the conventional, 2-grid IEC configuration. 
As illustrated in the first chapter, at the pressures which yield the highest neutron rates in 
the literature (4.0 microns) [18], the mean free path of ions in the system is less than the 
typical device diameter. The majority of these reactions were shown to be beam- 
background reactions [357]. As the pressure is reduced, the reaction rate peaks at 
approximately 0. 1 microns and then is also reduced, but the ion confinement time is 
increased. In the limit of a perfect vacuum, the only means of ion loss from the IEC 
system are: 

1) Ion neutralization on a grid wire via unconfined/chaotic trajectory or ion-ion 
collisional scattering 

2) Ion neutralization on the anode via ion-ion collisional energy spread 

3) Fusion 

Hirsch estimated the average number of passes through the center of a 2-grid device in 
this hard vacuum regime with a simple grid transparency argument [21]. 

8 = — Equation 3.20 

1 - v 

In this simple model, vis the grid transparency (typically 0.9 to 0.98) and 5 is the 
expected number of passes through the center of the device. The sole loss mechanism 
implicit in this model is number 1) in previous list, i.e. loss to the grid wires. 
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OOPIC was used to check the validity of this assumption. A simple 2-grid model was 
developed which did not contain any of the realistic field asymmetries which had already 
been shown to be capable of reducing the ion confinement time to 1-5 passes. Due to the 
finite grid size the grid wires are larger than they typically are experimentally which 
limits the modeled transparency. In this model the grid transparency is 88% which 
should yield 4 passes through the center on average according to the Hirsch model. The 
OOPIC model presented here assumes that the ions are Ar +1 and the cathode potential is - 
7500 V resulting in an ion bounce time of approximately 8 ps or a single pass time of 4 
ps. The average confinement time is therefore expected to be 16 ps. 



Figure 34: Low pressure, 2-grid OOPIC simulation, ion injector on 


Figure 34 shows a 2-grid OOPIC model with no electric field asymmetries and no neutral 
background gas. When the simulated ion injector is turned off, the number of ions in the 
system will decrease as ions are lost. The 1/e time of the decay curve can be compared to 
the expected confinement time based on the Hirsch formula. 
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Figure 35: Time evolution of Argon +1 ions in the simplest 2-grid OOPIC model (no background gas) 

Figure 35 shows the time evolution of ions in this model under these idealized conditions. 
Not surprisingly, when injection is turned off at 200 ps (top two images is figure 35), 
there is a rapid loss of ions on the order of the Hirsch confinement time. The top two 
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images in figure 35 show the spatial location of the ions on the left and the spatial density 
of the ions on the right. The rapid loss of ions appears to be localized in space - it is 
only the ions outside of what may be considered a “grid acceptance window” which are 
lost on this rapid 16 ps timescale. 

The middle two images in figure 35 show the ions in the system 230 ps after termination 
of injection. The image on the right is the same spatial representation of the location of 
the ions, while the image on the left shows the time history of the total number of ions in 
the system. It is clear that the rate of decay of ions in the trap has changed in such a way 
that it is no longer proportional to the number of ions in the trap as would be expected if 
confinement were either limited by statistical grid wire interception or collisions with 
background gas. In fact, the 1/e decay time of those ions inside the grid acceptance 
window is approximately 40ms — orders of magnitude longer than the 16 ps predicted 
from the transparency model. Even though the local field around the cathode is 
nominally defocusing, the overall ion confinement is good enough in this idealized 
simulation that the two-stream instability has time to develop and bunch the ions as seen 
in the bottom two images of the figure 35 at a time of approximately 1ms after 
termination of injection. This bunching was first identified by McGuire [357] in 
computational modeling of IEC devices and is a saturated mode of the well-known 
streaming instability in this type of device. 

When asymmetries in the field structure are modeled by introducing the realistic effects 
of a stalk or feed-through, however, the Hirsch model does much better at predicting ion 
confinement in a 2-grid IEC device. 
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Figure 36: 2-Grid OOPIC model showing ions impacting cathode after 1/2 to 5 passes 
Figure 36 shows how the field asymmetries caused by the cathode feed stalk result in 
very short ion lifetimes of 1-5 passes through the inner cathode grid — dependent upon 
the location of the ion injector(s). 

These simple simulations suggest that the “transparency” model developed by Hirsch 
may not be a good means of estimating the confinement time in a 2-grid IEC device if the 
device is carefully constructed to minimize the effects of asymmetric background fields, 
but it should be reasonable for most 2-grid devices with realistic field asymmetries due to 
feed through stalks. 

3.4.2. Multi-Grid Modeling and Confinement Estimations 

One very effective way of shielding the effects of feed-through asymmetry and 
improving ion confinement is by carefully tailoring the field structure with multiple 
independently biased grids. The efficacy of this approach was clearly illustrated using 
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OOPIC in chapter 1 figures 3 and 4. McGuire has conducted extensive investigations 
into this technique using OOPIC [357] . His original code was used as the basis for all 
OOPIC modeling of multi-grid IEC devices contained in this thesis. 

If we assume that the multiple grid configurations are capable achieving the highly 
confined ion trajectories observed in the simulations, the primary limiter on ion 
confinement is scattering due to collisions. Scattering with background gas and other 
ions and thermalization become the primary loss mechanisms. Therefore it should 
become possible in UHV conditions to push the theoretical limits proposed by Rider 
[175], 

In order to achieve these theoretical improvements, an experiment was designed with the 
goal of demonstrating the potential to significantly improve confinement time using 
multiple grids. This experiment is operated under high vacuum conditions; lower than 
the pressures at which most IEC experiments are run, but above the UHV pressures that 
would be required to push the limits of thermalization. At these pressures, collisions with 
the background gas are the biggest limiter to confinement. OOPIC was used to estimate 
the effects of background gas on ion confinement times in the multi-grid experiments. 
Figure 37 shows a typical, multi-grid simulation geometry in OOPIC. 
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number(t) 


x-y phase space for ions 



Figure 37: Typical simulated multi-grid geometry, x and y dimensions are meters, grid locations and 
ion injector geometry are based upon the physical experiment 

In order to facilitate evaluation of the 1/e time, the ion signal is norm a li z ed to the level at 
the time when the injection is turned off. This practice will be used throughout the thesis. 
The complete, un-normalized curves typically resemble the curve shown in figure 38: 


number(t) 



Time 


Figure 38: Typical plot of number of macroparticles in an OOPIC simulation 


The standard practice that was adopted for both computational experiment and physical 
experiment was to inject ions for 1ms then turn off the injector and observe only the 
decay portion of the curve. While the maximum number of ions in the system (and hence 
the signal strength) can vary with the strength of the source, the decay curve will only be 
a function of the sink. Comparing the 1/e times of the normalized decay curves will 
therefore give the best indication of actual confinement time of ions in the device. 

Figure 39 shows a map of the potential structure in a typical multi-grid OOPIC 
simulation which illustrates the ability of the multiple grids to “shield out” the effects of 
the feed-through stalks. 



Figure 39: OOPIC multi-grid potential model for a scenario with good ion confinement 

Not all multi-grid scenarios showed improved confinement. In fact the ion confinement 
was expected to be quite sensitive to the particular grid potentials based upon OOPIC 
modeling. 
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Figure 40: Symmetric multi-grid field with poorly confined ions 
In the figure 40 above, the field structure is symmetrical, but the ions are not well 
confined due to a field structure that prematurely focuses the ion beam. The beam then 
freely expands in the core region where there is little confining field curvature. On the 
opposite side of the device, the over-focused field again skews the trajectory of the beam 
which results in a wider beam envelope. Ions are quickly lost to the grid wires, and the 
overall confinement is very poor. It is clear from this one example that the specific 
potential structure imposed by the grids can be quite important in determining the ion 
confinement time even in the absence of field asymmetries. 

When a small asymmetry in one of the grid wires is introduced (a one grid space 
perturbation which corresponds to a 2.5% change in position of one grid wire), the poor 
confinement is even more pronounced. Figure 41 shows the effect of an unperturbed 
(left) and perturbed (right) wire in the multi-grid simulation. The perturbation introduces 
field asymmetries which tend to skew the ion beam 


82 



Figure 41: Comparison of perfectly symmetric field (left) with a 2.5% location perturbation on one 
wire (right). The perturbed wire is the second grid from the center in the upper left hand quandrant. 


If the imposed potential structure results in a well confined ion beam (figure 42), 
however, there is much less sensitivity in the overall shape of the beam envelope to 
geometric perturbation in either the grid wires or the ion source. 



Figure 42: Good multi-grid confinement with unperturbed grids and concentrated ion source (left), 
and highly perturbed grids (5% on 2 nd grid and 2.5% on 3 rd grid) with diffuse ion source (right), 
same grid potentials 


The geometry on the right hand side of figure 42 above is representative of the worst case 
scenario for the experiment as designed at an achievable pressure of 8e-5 torr. For this 
model with Helium gas, He +1 ions and a cathode potential of -5000 Y, the predicted 1/e 
confinement time is 8 ps or approximately 2 passes as shown in figure 43. 
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Normalized Ion Content, He, 8e-5torr, worstcase 



Figure 43: Normalized plot of the decay of the number of ions in the "worst case" OOPIC simulation 
for 8e-5 mbar He. 

Once the experiment was completed and data with various diagnostics was taken, it 
became clear that for Helium, the pressure range with the highest S/N was le-4 mbar or 
8e-5 torn A comparison of best expected confinement time at this pressure was then 
conducted to evaluate how close to the best simulated confinement we had achieved in 
the actual experiment. The worst case scenario is presented in figure 43 above. 

Although the worst case 1/e confinement time of 8ps is quite bad, it is still better than the 
perfect vacuum 2-grid simulation (5.4 ps) with the stalk assymmetry. 

The decay for the best case scenario (left side of figure 42) is shown in figure 44 below. 
The grid potentials from anode to cathode for both of these runs are [0, 0, -1000, -1500, - 
5000], It is worth noting that the geometry of the ion source combined with large 
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(though not unreasonable) asymmetries in the grid wires have the potential to alter the 1/e 
confinement time of the system by a factor of 5. This change is purely geometrical as a 
result of the induced field asymmetries, i.e. the pressure, the gas, and the grid potentials 
are all the same in these two scenarios. 


Normalized Helium Ion Population Decay Curve, 8e-5torr 
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Figure 44: Normalized decay curve of best case, Helium at 8e-5 torr 
These results can also be compared with estimations of the mean free path for collision 
and charge exchange at this mean ion energy and device pressure. 


3 1 _ kT 


Equation 3.21 


Assuming a room temperature background gas at le-4 mbar, the mean free path for 
charge exchange at this pressure and average energy is only 4. 1 m corresponding to 8 
passes or 32 ps based on the charge exchange cross-sections given by Hegerberg [352], 
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The OOPIC simulation in conjunction with these simple calculations suggests that the 
dominant loss of ions at this pressure is due to charge exchange collisions with the 
background gas. While all charge exchange reactions will result in some loss of energy 
(via the fast neutral), not all charge exchange reactions result in direct loss of an ion from 
the system (especially those ions that undergo charge exchange near the edge of the 
potential well), it is therefore possible that confinement times longer than the mean time 
until charge exchange could be expected. 

The effects of changing the background gas pressure were simulated. For all of these 
multi-grid runs, the same confining potentials were applied to the grids [0, 0, -1000, - 
1500, -5000] volts, and the corresponding grid radii are [.25, .17, .10, .075, .05] meters. 
As expected, the 1/e confinement time is proportional to 1/p for multi-grid experiments. 
1/e confinement time for 2-grid experiments shows a general insensitivity to pressure in 
the regimes explored. Figure 45 clearly shows the expected difference between 2-grid 
and multi-grid confinement times based on OOPIC modeling. 
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1/e Confinement Time v. 1/Pressure 



Figure 45: OOPIC simulated effect of pressure on 1/e confinement time for 2-grid and multi-grid 
IKC systems 

It is worth noting that the simulated 2-grid confinement times are implicitly worst case 
scenarios. The 2-dimensional simulation results in a loss of all particles that pass by the 
stalk. In the real system, not all ions with “low” trajectories will impact the stalk on then- 
first pass. Unfortunately full 3D PIC modeling of this system was not available for these 
investigations. 

Figure 46 clearly shows the expected difference in the ion trajectories in the multi-grid 
(left) and 2-grid (right) geometries. 
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Figure 46: Multi-grid (left) and 2-grid (right) OOPIC simulations 


Simulated runs at zero background pressure were also conducted for the multi-grid and 2- 
grid scenarios above. At zero pressure, the 2-grid 1/e confinement time was 5.4 ps (still 
limited to -lpass), and the multi-grid 1/e confinement time was 2640 ps (2.6 ms). In the 
zero pressure limit, the loss rate of ions from the multi-grid system should go as n not 
simply n, so it may be desirable to characterize confinement times with a different metric 
than the 1/e time used in this work. This is due to ion-ion collisional effects as opposed 
to ion-background collisions which dominate in the current experiment. A difference in 
the shape of the decay curve was noted in the simulation (figure 47), which would 
suggest this relationship. Due to the very long confinement times and hence the very 


long simulation run times, this relationship was not thoroughly explored. 


number(t) 



Time 


Figure 47: Non-exponential decay shape of multi-grid in perfect vacuum 
In order to achieve this different confinement regime in an experiment of this size where 
the density of ions is greater than or on par with the background gas density assuming a 
comparable strength ion source, base pressures on the order of 10' n mbar would need to 
be achieved. This n confinement relationship cannot be explored in experiments at this 
time due to the limitations of the existing vacuum equipment and ion sources. The multi- 
grid UHV regime exploration is left to future students. 


In addition to the pressure sweeps shown in figure 45, a simulated voltage sweep of the 
3 rd grid was also conducted in order to evaluate the sensitivity of the 1/e confinement 
time to the potential structure. The potential on the third grid was varied from 0 (the 
potential of the 4 th grid) to -1500 V (the potential of the second grid). 
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OOPIC Simulations 



-1600 -1400 -1200 -1000 -800 -600 -400 -200 0 


3rd Grid Potential (Volts) 

Figure 48: Predicted sensitivity of 1/e confinement time to 3rd grid potential (le-4 mbar) 

It can be seen in figure 48 that the particular grid potential can have a significant impact 
on the confinement time of ions in the system. The purpose of this sweep is to compare 
the OOPIC modeling to the experiment so the effectiveness of the 2-D PIC model can be 
determined for predicting confinement behavior. Experimental results will be presented 
in the next chapter. 

3.4.3. Summary of Expectations Based on Modeling 

Due to the presence of field asymmetries, all 2-grid experiments are expected to have ion 
confinement times roughly predictable by the Hirsh transparency model. The accuracy of 
that model is questionable, but in general, 2-grid experiments are not expected to show a 
significant dependence of confinement time on pressure because confinement will be 
limited by skewed trajectories resulting from asymmetric potential fields. At a -5000 V 
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cathode potential with Helium ions, the expected 1/e confinement time is 28 ps (7 passes) 
while Argon ions should be confined for approximately 84 ps. 

Multi-grid experiments with proper tuning of the grid potentials, however, are expected 
to show a pressure dependant confinement time because of their ability to shield out field 
asymmetries. It is expected that the multi-grid design will result in a significant 
improvement in 1/e confinement time over the conventional 2-grid IEC configuration. 
OOPIC modeling has suggested that at a pressure of 8e-5 torr (He), with a -5000Y 
cathode, the 1/e time of ions in the system could be 43 ps - significantly better than the 
best 2-grid estimate at these potentials. 

If improved confinement is achieved, there may be signs of the saturated two-stream 
instability as has been seen in the idealized OOPIC simulations. 
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4. Experimental Results and Discussion 


This section gives an overview of the data collected on the MIT SSL IEC experiment. 
Initial testing was done with only the 10 cm diameter cathode grid and the 50 cm 
diameter anode grid. The three intermediate grids (15, 20, and 35 cm diameters) were 
then added for multi-grid testing. Due to changes in the ion source and confinement 
diagnostic during the multi-grid testing, a final round of 2-grid testing was done to have a 
more direct comparison. 

For both configurations, testing was conducted under low and high vacuum conditions. 
The low vacuum testing (higher pressure) was conducted in air while the high vacuum 
data used a backfill of Helium or Argon gas depending upon the type of ion desired. All 
low vacuum data acquisition was done with a 6 megapixel digital camera. High vacuum 
data acquisition was done with a 2 GHz Tektronix TDS 2014 Oscilloscope. 


4. 1. The 2-Grid Configuration 

Data was first collected in a traditional 2-grid device. The simple anode-cathode 
structure (as shown in figure 1) is the most abundant IEC Fusion reactor design and has 
been extensively studied by Farnsworth, Hirsch, Miley and many others. The data taken 
in this configuration would serve as a baseline for purposes of comparison. Figure 49 
shows a picture of the 2-grid device in the vacuum chamber. (Note: the cylindrical 
structure in the port on the left of the image is an old ion source that was not used). 
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Figure 49: 2-Grid assembly in the vacuum chamber 

4.1.1. Low Vacuum Operation 

Early experiments were done in the low vacuum regime in order to satisfy our curiosity 
about the optical discharge that would appear. 


The background pressures during these discharges are sufficiently high that ion-electron 
cascading results in a rapid change in discharge impedance. As more ions are generated, 
the impedance decreases and the current rapidly begins to grow as secondary electrons 
liberated from the grid wires by ion impact cause more ionization of the background gas. 
The power and current limits of the high-voltage power supplies connected to the grid 
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imposed by the built-in protection circuitry result in a highly irregular pulsing of the 
voltage and current on a time scale of a few milliseconds. Because these voltage and 
current fluctuations are driven by the power supply protection circuitry instead of the 
actual physics of the plasma discharge, the shape of the voltage and current pulses are not 
worthy of study. A consequence of this cascading phenomenon is that the most 
interesting “real” data are the measurements of the optical discharges in the various 
pressure regimes as captured by digital camera. 



Figure 50: Jet-mode in 2-grid configuration 

It was found that all of the reported modes of operation of high pressure IEC devices 
could be successfully recreated. Namely, there was evidence of both the jet-mode plasma 
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configuration as shown in figure 50, and the star-mode plasma as shown in figure 51. 
These modes have been previously reported by Miley et. al.[155]. 



Figure 51: Star mode (1.7e-3 mbar) in 2-grid configuration 


4.1.2. HV Operation 

The low vacuum operations were interesting from an aesthetic point of view and to 
confirm results reported by others, but the primary interest was increasing the ion 
confinement time which could only be done in high vacuum. For most of these 
experiments, a base pressure of approximately le-6 mbar was reached, and the chamber 
was then back-filled with the gas of choice for ionization. 
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4.I.2.I. Measurements of current to the grid wires 

Early experimentation was done in the two grid configuration at low voltage (-400Y 
cathode) to evaluate the possibility of using the measurement of the currents to the grid 
wires as a measurement of ion confinement time. A simple voltage divider circuit was 
used to measure the small fluctuations in grid potential that indicate the ion flux to the 
grid wires. The detection limits of the circuit were tested by simulating a high impedance 
pure current source using a high voltage square wave signal generator in series with a 
lOMOhm resistor attached to the cathode grid. The RC decay of the detection circuit can 
be seen in figure 52 to be approximately 3 microseconds, much shorter than the shortest 
measurable confinement time in the system (Helium ions with an obstruction). 


He Decay and Detection Limits 



Time (s) 


Simulated Current Source Test 36: He Large Foil Test 35: He, No Obs 


Figure 52: Detection limits of sensing circuit due to RC smoothing 
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Helium Decay Current, Varying Voltages 



Test 35: -400V (-344 to -358V), 1 ,4e-4mbar He Test 37: -200V (-1 82 to -186V), 1 ,3e-4 mbar He ] 

Test 38: -1 00V (-86 to -90), 1 .4e-4mbar He | 


Figure 53: Effect of varying cathode voltages on Helium ion confinement time 
Figure 53 shows that when the cathode potential was varied, there was an effect on the 
confinement time. In fact these data confirmed that reducing the cathode potential has 
the expected effect of reducing the 1/e confinement time varying roughly with the inverse 
of the square root of the absolute magnitude of the cathode potential. This supported the 
argument that the ion confinement was limited to a certain number of passes through the 
system. Due to the experimental limitations on this detection technique discussed in 
chapter 2, this technique will not be discussed further. 

4.I.2.2. Measurements with the capacitive probe 

While measurements of current to the grid wires gave a good initial grounding to these 
investigations, the capacitive probe allowed much lower cathode potentials to be 
achieved. The advantages of this measurement technique were expounded upon in 
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chapter 2. This section presents the data collected in the 2-grid configuration with the 
capacitive probe under a wide range of pressure conditions. 


2-Grid, Helium, Pressure Effect 



Figure 54: Effect of pressure on 2-grid confinement. Helium 

Figure 54 shows the 2-grid helium confinement time plotted against the inverse of 
pressure. These data clearly indicate that pressure does effect confinement time in 2-grid 
IEC devices. They also show what appears to be two distinct confinement regimes, a 
high pressure regime and a low pressure regime. The dividing inverse pressure is around 
5000 mbar' 1 or 2e-4 mbar. At pressures higher than this dividing pressure, the 
confinement appears reasonably linear with 1/P as shown in figure 55. The slope of the 
linear regression is similar although slightly lower than the slope of the predicted multi- 
grid confinement, and the regression fit is quite good. 
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At pressures lower than this dividing pressure (figure 56), there is no clear relationship 
between confinement time and pressure and the regression fit is very poor but it shows 


that the confinement time is does not increase as the pressure is further reduced. 


2-Grid, Helium. High Pressure Regime 





Figure 55: 2-Grid confinement, high pressure regime, Helium 
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2-Grid, Helium, Low Pressure Regime 



Figure 56: 2-Grid confinement, low pressure regime, Helium 
This insensitivity of confinement time to pressure is similar to what was predicted with 
the OOPIC modeling, although the OOPIC modeling showed a much lower confinement 
time than is measured in experiment. Possible reasons for this discrepancy include the 2- 
dimensional nature of the OOPIC model and the lower simulated transparency of the grid 
and stalks. The ion slug at injection termination also artificially increases all of these 
numbers equally by approximately 7-9 ps for helium. 

The 2-grid Argon data (figure 57) show a similar 2-regime effect, although the critical 
pressure for Argon appears to be around le-4 mbar. The ion slug artificially increases 
these measured confinement times by approximately 10-12 ps for argon. 
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2-Grid Argon 



Figure 57: 2-Grid confinement time pressure sensitivity, Argon 

The Argon data indicate that there may be a very slight effect of pressure in the very low 
pressure regimes, but it is clearly not as strong as it is at high pressure; the slopes differ 
by a factor of 22 between the two regimes as shown in figures 58 and 59. 
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2Grid Argon, High Pressure Regime 



Figure 58: 2-Grid high pressure regime, Argon 


2Grid Argon, Low Pressure Regime 





Figure 59: 2-Grid low pressure regime, Argon 
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The experimental data collected in the 2-grid configuration with both Helium and Argon 
gas indicate that the 2D OOPIC model does not yield an accurate prediction of ion 
confinement time in a 2-grid IEC device. 


A simple analysis shows that the relative 1/e confinement time trends of Helium and 
Argon in their respective high pressure regimes are consistent with a model in which the 
primary loss mechanism is due to charge exchange with background gas. 


1/ e _ConfinementTime oc 

IonLossRate 

IonLossRate oc n i n h a cx v i 

1/ e _ConfinementTime oc — - — 

n b a cx V i 

Therefore, at the same pressure (background density), 

1 /e_Ar <J cxHe v He 

1 !e_He cr cxAr v Ar 


Equation 4.1 
Equation 4.2 
Equation 4.3 


Equation 4.4 


The charge exchange cross-section of Argon at the energies of interest is approximately 
2.6 times the cross-section of Helium [372], but the 10: 1 mass ratio for the two species 
results in a theoretical 1/e confinement time ratio of 1.21 for Ar/He. The data presented 
in this section shows a 1/e confinement time ratio of 1.38 at a pressure of 2.2e-4 mbar 
(the low end of the Helium high pressure regime and the high end of the Argon high 
pressure regime data). This 14% deviation can be explained by a combination of the 
inaccuracy of the data collected here due to the “slug effect” and imprecision in the 
literature on charge exchange cross-sections (~5-10% spread in the data), the assumptions 
of the energy at which most charge exchange occurs that results in ion loss, and finally 
(and perhaps most importantly) the influence of the presence of the grids with a finite 
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transparency. Since Argon is more massive, the ions will have ~l/3 rd the “passes” 
through the system in a given amount of time as Helium. Helium consequently has many 
more opportunities to hit a grid wire. This could help explain the larger measured spread 
of confinement time than predicted by the pure charge exchange model. It could also 
help explain why the dividing pressure for Helium is much higher than the dividing 
pressure for Argon. 

The flattening of the trends at lower pressure indicates that the ion lifetime is limited 
substantially by a mechanism other than charge exchange such as direct grid impact of 
primary ions. This would be consistent with the flat confinement time predictions from 
OOPIC or the Hirsch model - a statistical model in which on average ions are lost by 
cathode impact after a certain number of passes through the system. In the zero pressure 
limit, one would expect ratio of 1/e confinement times of Argon/Helium to approach the 
square root of the mass ratio (3.16). These data indicate a ratio of 1.5 to 3 in confinement 
time at the lowest measured pressures (the error bars are large at the lowest pressures). 
These data also suggest that the number of passes at the lowest pressures is limited to 
approximately 10 ion-core transits, quite a few more than predicted by the Hirsh 
transparency model. 
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4.2. Multi-grid Configuration 


The core hypothesis of this thesis is that ion confinement times can be increased by the 
use of multiple (>2) independently biased grids. This section reports on the first known 
multi-grid IEC experiment. 

4.2.1 . Low Vacuum Operation 

Again, the device was first tested at low vacuum in order to see the optical discharge. 
Interesting plasmas were generated. They are presented here for purposes of 
completeness and qualitative comparison with the 2-grid pictures. 



Figure 60: Multi-grid discharge 2e-2 mbar, air 

At very high pressures, the discharge occurred at the feed through. As the pressure is 
dropped, the discharge appeared at the polar regions as in figure 60. As the pressure 
continued to drop, distinct jets became visible (figures 61 and 62) and gradually grew 
more pronounced at the equator. 
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Figure 61: Many jets (~1.5e-2 mbar) 



Figure 62: Multi-grid air discharge le-2 mbar 
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Figure 63: Multi-grid jet-mode (7e-3 mbar) 





Lower pressures resulted in a single primary jet that would move between equatorial 
longitude lines as shown in figures 63 and 64. These jets were a more bluish color than 
the higher pressure pink jets. Eventually, as the pressure was lowered still more, the 
optical discharge would reduce in intensity and the star-mode (figure 65) with its 
characteristic multiple beams and small core focus was developed (although at a slightly 
higher pressure than in the 2-grid runs). The discharge also terminated at a higher 
pressure (2.7e-3 mbar) than it had in the 2-grid runs (1.5e-3). 



Figure 65: Multi-grid, star-mode 3e-3 mbar, air 


It is worth noting that during these high pressure discharges, the potentials on the grids 
are not held constant. The HV power supplies on the grids have current limiting circuitry 
which limits the potential and the discharge on a 1ms time-scale. The photon capture 
time for the CCD on the camera used to take these pictures is up to 0.6 seconds, so these 
rapid pulses are not seen. In fact, the “double -jet” image in figure 64 is somewhat 
deceiving; the jet moved between those two locations w hi le the electronic shutter of the 
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camera was open. There was only one jet at any given time in that regime. In the star- 
mode regime, the jets appear to be simultaneous to the resolution of the naked eye 
(<~0.1s). 

4.2.2. High Vacuum Operation 

Base pressures of 5e-7 mbar to 1.5e-6 mbar were typically reached for the multi-grid high 
vacuum experiments. 

4.2.2.I. Effect of Pressure on 1/e Confinement Time 

In Chapter 2 interactions of ions with the background gas was identified as the primary 
limiter on ion confinement time in multi-grid systems with good confining potential 
structures. It was therefore desirable to experimentally evaluate the effect of pressure on 
1/e ion confinement time. For the following data, a potential structure with good 
confinement was imposed by setting the grid potentials at [-5000, -1500, -1000, 0, 0]. 

The background pressure was varied using a fill gas of Helium and Argon. Since 
confinement time is expected to increase linearly when plotted against the inverse of 
pressure, that is the format that has been adopted herein. 
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Effect of Pressure on Confinement 



Figure 66: 1/e Confinement time v. 1/Pressure, Helium 


Pressure Effect: Comparison of Experiment to OOPIC 
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Figure 66 shows the Helium confinement time plotted against the inverse of pressure. 
Figure 67 compares those data to the OOPIC predictions for the same potential structure 
over a similar range of background pressures. It is interesting to note that the slopes of 
both the physical and computational experiments are the same, but the data from the 
physical experiments show a longer confinement time than predicted by the OOPIC 
modeling. This constant 20 ps offset can be explained by a combination of three 
differences between the OOPIC model and the experiment. 

1) The modeled transparency of the grid wires in the OOPIC simulation is only 87% 
due to computational grid limitations but the actual, physical transparency is 93%. 

2) The leakage of the slug of ions from the ion source outside the anode is not 
modeled in the OOPIC simulations. 

3) The OOPIC model measures the total number of particles in the system which 
starts decaying immediately when the ion injection is stopped whereas the 
physical experiment detects the ion density near the probe (opposite the ion 
source), so there is a delay corresponding to the ion transit time across the system 
(~4 ps for He +1 and 12 ps for Ar +1 for a -5000 Y cathode) 

Of the unaccounted for 20 ps, 4 ps is clearly explained by the delay due to the location of 
the probe, the remaining 16 ps can easily be accounted for by a combination of the 
transparency factor and the ion slug. It is not clear how much of the remaining 
discrepancy is due to which factor since the velocity and birth location of the average ion 
in the slug can vary significantly, and the magnitude of the field penetration from the 
cathode cannot be accurately predicted with the tools at hand, but the data consistently 
indicate that there is a peak in the ion signal 16-17ps (figure 68) after termination of 
injection in these multi-grid tests. 
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Helium, 1.0e-4 i 



Figure 68: Multi-grid Helium ion decay curve, le-4 mbar 
Compensating for this 16-17 ps lag in combination with the 4 ps ion transit time would 

result in an exact overlap of the two curves in figure 67 to within lps over the entire 
range. This would result in the linear data fit of figure 67 passing through the 0,0 point 
as would be expected. 

Figure 69 shows a similar peak in the ion slug for Argon ions. It is interesting to note 
that the initial magnitude of the “slug effect” is comparable between Argon and Helium, 
but the dampening of the slug signal is faster at the same background pressure with 
Argon than with Helium (likely due to Argon’s higher charge exchange cross-section). 
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Argon, 1.0e-4mbar 



Time (s) 


Figure 69: Multi-grid Argon ion decay curve, 1.0e-4 mbar 


The OOPIC model in combination with these simple data analyses is thereby shown to 
give a reasonable approximation of the effect of pressure on the confinement time of ions 
in a multi-grid IEC device. Although exact numerical prediction was not achieved due to 
the transit time and the slug effect, the slope of the variation of the 1/e confinement time 
with the inverse of the pressure was shown to be within 5% of the OOPIC prediction for 
Helium. The Argon data is even more interesting. 
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Multi-grid Argon Confinement Pressure Sensitivity 



Figure 70: Multi-grid Argon Confinement Pressure Sensitivity 

It can be seen in the figure 70 above that the data for argon do not provide as good a 
linear fit as the Helium data. The slope of the fit is also nearly half of the slope of the 
Helium data. Due to Argon’s much larger electron-impact ionization cross-section (a 
factor of ~ 8), an acceptable signal to noise ratio is maintained at much lower pressures 
with Argon than with Helium. Note the inverse pressure range for Helium is 1 GOO- 
13, 000 mbar' 1 in figure 67 while the range for Argon in figure 70 is 5000-25,000 mbar' 1 . 
It was desirable to expand this range as much as possible given the available equipment. 


An interesting phenomenon began to occur around le-4 mbar in the experiments with 
Argon: the slug of ions that is emitted when the ion source is turned off results in an 
underdamped ringing at pressures lower than around le-4 (figure 72). At pressures 
higher than le-4, this ringing is damped within a few cycles (figure 71). 
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Argon 1.33e-4mbar 



Time (s) 


Figure 71: Higher pressure regime decay example 
Argon 6.8e-5 mbar 



Figure 72: Lower pressure regime decay example 
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While this ringing signal appears similar to the ringing one would expect from the two- 
stream instability, it is a direct result of the ion slug emitted at injector shut-down, and 
should not be interpreted as direct evidence of the instability. It appears, however, that 
when this “ringing” is present, the ion confinement does not follow the same trend as is 
established at higher pressures when the “ringing” is damped out within a few cycles. 
Indeed, if the data in figure 70 is split into two separate regimes, one at higher pressures 
where the ringing is damped, and one at lower pressures where the ringing is highly 
underdamped, linear fits in both regimes are much better. 

ArgonPressureEffect_highpressregime 



Figure 73: High pressure confinement regime, Argon 
In the high pressure regime (figure 73), the slope of the linear fit is similar to the Helium 
data to within 12%. The linear fit as measured by R 2 is also quite good — comparable to 
the fit of the Helium data in the similar pressure regime. 
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Low Pressure Regime Only 



Figure 74: Low pressure confinement regime, Argon 

The data in figure 74 suggest that the presence of an underdamped, “ringing” bunch of 
ions could significantly reduce the confinement time based upon the comparison of the 
slopes of the linear reductions in the two regimes. The error bars in the lower pressure 
regime are much larger due to the difficulty in ascertaining the exact 1/e time in the 
presence of a proportionally large amplitude oscillation. 

4.2.2.2. Intermediate Grid Voltage Sweeps 

It was desirable to evaluate the sensitivity of the ion confinement to the potential 
structure. For all of the tests reported in this section, the 5cm radius (1 st ) grid is held at a 
potential of -5000 V, and the 25cm radius anode grid is held at ground. 
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The baseline grid potentials are [-5000, -1500, -1000, 0, 0]. The second, third, and fourth 
grid potentials were varied individually from this baseline condition in the experiment to 
evaluate the effect on ion confinement time. 


3rd Grid Voltage Sweep, P = 1e-4 mbar, He 



-1600 -1400 -1200 -1000 -800 -600 -400 -200 0 


Grid Voltage 

■ Tests 178-201 — OOPIC Simulations Poly. (Tests 178-201) 

Figure 75: Voltage sweep comparison to OOPIC simulation 
The 3 rd grid or middle grid (20 cm diameter) was expected to have less of an effect on the 
confinement of ions in a multi-grid system than either the 2 nd or the 4 th grids based upon 
extensive computational and physical testing. It was therefore desired to use a sweep of 
the potential of this grid to evaluate the ability of the 2D particle-in-cell code, OOPIC, to 
predict the confinement of ions in a 3D physical experiment. As shown in figure 75, The 
ability of OOPIC to predict the effect of the variation of grid potential on the ion 
confinement time was shown to be quite poor based on these data. OOPIC predicted that 
in the range of -600 to -1200 volts, there would be a substantial change in the 
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confinement time. No such significant change was visible in the data until the 0 to -400 
volt range was examined. It is clear based upon these data that the 2D OOPIC code does 
not give an accurate prediction of the effect of the 3 rd grid potential on the confinement 
properties of the multi-grid IEC system in the experiment. 

The only conclusion that can be reasonably drawn from these data is that OOPIC 
consistently seems to under predict the real, measured confinement times, and as the 3 rd 
grid potential approaches the potential of the next grid in, the ability to predict the 
confinement seems to improve. Some possible reasons for this under prediction were 
expounded upon earlier in this section. One possible reason that was not explained 
before is the slight difference in the field structures due to the 2D Poisson solver versus 
the real, 3D scenario. This slight difference can affect the confining properties of the 
potential structure. The additional degree of freedom in experiment also can allow 
particles to move around things like feed-through stalks in the real world whereas in the 
2D simulation, all particles going “past” the stalk, hit the stalk. It would be quite 
informative to compare the experiment to a full 3D PIC code. 

Due to its inability to accurately predict the sensitivity of ion confinement to intermediate 
grid potentials, OOPIC modeling is not presented with the following grid sweep data. 
Figures 76 and 77 show the effect of 2 nd grid voltage (15 cm diameter) on the 
confinement time for Helium ions and Argon ions respectively, Figures 78 and 79 show 
the effect of 3 rd grid voltage (20 cm diameter) on the confinement time of Helium and 
Argon ions respectively. And figures 80 and 8 1 show the effect of 4 th grid voltage (35 
cm diameter) on the confinement time of Helium ions and Argon ions respectively. 
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Voltage Sweep of 2nd Grid, 9.4e-5 mbar, Helium 
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Figure 76: 2nd Grid voltage sweep, Helium 

Voltage Sweep 2nd Grid Argon, 9.2e-5mbar 
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Figure 77: 2nd Grid voltage sweep, Argon 
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Voltage Sweep of 3rd Grid, 1e-4 mbar Helium 
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Figure 78: 3rd grid voltage sweep, Helium 
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Figure 79: 3rd grid voltage sweep, Argon 
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Voltage Sweep 4th Grid, 7.7e-5 mbar Helium 
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Figure 80: 4th grid sweep, Helium 
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Figure 81: 4th grid sweep, Argon 
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The data on the previous pages clearly show that the sensitivity of ion confinement to 
grid potential is not strongly dependent on the type of ion. While argon confinement 
ti mes were all generally longer than Helium confinement times as expected, the 
sensitivities to grid potentials were consistent between the two ion species. It can also be 
concluded from these data that ion confinement is most sensitive to the outermost grid 
potentials i.e. those near the ion “turn-around” potential. The fourth grid clearly exhibits 
the most sensitivity to grid potential. 
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4.3. Summary of Data 


This section summarizes the most important results obtained from the experiment. Each 
data point represents an individual experiment. For brevity, the raw data is not included 
in this thesis except in a few instances for purposes of extreme examples at very low 
pressures. 

4.3.1. Comparison of 2-grid to multi-grid data 

The primary goal of this research was to show that multiple independently biased 
concentric spherical grids could be used to improve ion confinement times in IEC 
systems. The following pages present the comparison of all of the data collected that can 
be used as a basis for this comparison. These data were all collected with the capacitive 
probe detector under high vacuum conditions with both Helium and Argon ions in a 
background of neutral Helium and Argon respectively. The ion source and detector 
probe were identical for all of these experiments. The multi-grid data were collected with 
five concentric spherical grids of radii 5cm, 7.5cm, 10cm, 17cm, and 25cm, and grid 
potentials of -5000V, -1500V, -1000V, 0V and 0V respectively. For the 2-grid tests, the 
three intermediate grids were removed and the 5cm grid was held at -5000V while the 
25cm grid was held at ground. The probe potential was held at 0.3V above ground, and 
the ion source was turned on and off by rapidly (~2ps) changing the filament potential 
from ground to -150V (on) and back to ground. 
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Pressure Effect, Helium 



Figure 82: Helium, Multi-grid v. 2-grid confinement 


Figure 82 compares the 2-grid and multi-grid data for all of the experiments conducted 
with Helium ions and Helium background gas. At high pressures, the 2-grid device has 
better confinement than the multi-grid device. This is likely due to the short mean-free- 
path of Helium ions relative to the size of the device and the larger ion sink surface area 
present in the multi-grid tests. As the pressure is lowered, below le-4 mbar, the multi- 
grid confinement time is seen to exceed the 2-grid confinement time by a significant 
margin. These data strongly suggest that a multi-grid IEC operating at very low pressure 
could have a significantly longer confinement time than a 2-grid IEC device. 


125 



Multi-grid v. 2-Grid, Pressure Effect, Argon 



« Multi-Grid Data 
■ 2-Grid Data 

Linear (Multi-Grid Data) 

Linear (2-Grid Data) 


Figure 83: Argon, Multi-grid v. 2-grid confinement 


Figure 83 compares the 2-grid and multi-grid data for all of the experiments conducted 
with Argon ions in Argon background gas. Again, it can be inferred from these data that 
a multi- grid IEC operating at very low pressure (below le-4 mbar) will have a much 
longer confinement time than a 2-grid device. 


It is worth noting, that although there exists both a high pressure regime and a low 
pressure regime in these multi-grid tests with Argon, similar to the 2-grid experiments, 
the slope of the low pressure confinement is a factor of 16 higher in the multi-grid data 
than in the 2-grid data. These data clearly indicate that a significant improvement in the 
ion confinement can be achieved at low pressures with a multi-grid IEC experiment. 
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It is clear that improved confinement has been achieved in experiment and the trends 
indicate that as pressure is reduced, the confinement will improve faster with multiple 
grids than with the conventional 2-grid device. Lower pressures were also explored 
particularly with Argon where the higher ionization cross-section allowed a discemable 
signal down to pressures as low as 1.9e-6 mbar. 



Figure 84: Comparison of Argon ion signals at very low pressures 

At these very low pressures, the error bars on the 1/e confinement time are very 
significant, and although the S/N is reasonable, the interpretation of the signal is difficult 
in the multi-grid system due to the presence of the ion bunching as seen in figure 84 
above. If instead of trying to discern the 1/e confinement time, we look for the time at 
which the S/N ratio is approximately equal to 1 (the time at which ions can’t be detected), 
the above raw data clearly show the potential of multi-grid IEC to improve ion 
confinement in that a clear ion signal is visible 1ms after termination of injection in the 
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multi-grid case, while the ion signal can only be seen over the background noise for 0.2 
to 0.3 ms in the two grid case. It should be noted that the data in figure 84 is raw, un- 
normalized probe voltage data which is why the signal is seen to “decay” in a positive 
voltage direction. 

4.3.2. Detection of two-stream instability 

A secondary goal of this work was to look for evidence of the two- stream instability that 
McGuire had seen in computational experiments with good confinement [357]. The 
saturated mode of this instability was generally seen to be either a single or double 
bunching of the ions as they “ring” in the well of the IEC system. This growth can be 
seen in the raw data of figure 84 and more clearly in figure 85 below. 


Multi-grid Anjou. 2.0e-6 inbar 



Time (s) 

Figure 85: Evidence of two-stream instability at the bounce frequency, Argon 1.9e-6 mbar 
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Figure 85 shows very clear evidence of the presence of instability. The detection of this 
single bunched mode was not difficult to achieve with the capacitive probe due to its 
sensitivity to the local ion density at the reflection point. Figure 85 above clearly 
illustrates the growth in amplitude of the oscillating signal after the ion injection is 
terminated (time 0). While a slug of ions is expected after the ion source is shut down, 
this bunching would gradually spread out (decay) over time or stay at roughly the same 
amplitude in the worst case, unless there is a collective mode instability. The only 
mechanism for an actual amplitude increase after all ion sources have been terminated is 
a collective instability. No doubt this instability is excited by the perturbation caused by 
the slug of ions as the phase of the signal matches the slug quite well. 

The instability was only clearly observed with Argon ions. This is believed to be due to 
Argon’s much larger electron bombardment ionization cross section. While a distictive 
S/N with Helium was only possible down to pressures of 4e-5 mbar, Argon signals could 
be resolved at pressures as low as 1.9e-6 mbar (figure 85). The instability was seen in 
argon from this lowest pressure up to 3.4e-5 mbar when the ion source was operating. 

As can be seen from the preceding figure, the bunching is very long-lived. In the test 
presented above at 1.9e-6 mbar, the ion signal can still be seen clearly 1.25ms after the 
termination of injection - indicating much better confinement than any other reported 
test, as would be expected at this low pressure. The data show 52 full cycles of the ion 
bunch or 104 passes of ions through the core of the device over a period of 1200 
microseconds. The corresponding cycle time of 23 ps is quite close to the bounce time 
predicted in OOPIC (~24 ps) for Argon ions. It can be clearly inferred from these data 
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that this is an Argon ion bunch oscillating at the bounce frequency. It is worth noting 
that the number of passes indicated by these data suggests that other effects such as the 
previously discounted magnetic drift may begin to have a significant impact on the ion 
confinement time in this long-confinement time operating regime. 


Normalized Oscillation Envelope 
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Data Curve Fit 

Figure 86: Normalized oscillation envelope and exponential-sigmoid curve fit for the test plotted in 
figure 85 (Argon, 1.9e-6 mbar) 

If we attempt to extrapolate the 1/e confinement time from figure 85 based on the 
amplitude of the oscillation envelope in the signal (figure 86), it can be seen in that the 
1/e time is approximately 400 or 33 passes. It is interesting to note that if the trend in 
figure 74 were to be extrapolated to this low pressure, the predicted confinement time 
would be 875 ps, so it is clear that the actual confinement time does not follow that linear 
extrapolation under these circumstances. It may be that the presence of the instability 
causes a more rapid loss of ions than would otherwise be expected. 
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The curve fit shown in figure 86 is a decaying exponential multiplied by a sigmoid 
function. The decaying exponential represents the continuous loss of ions from the 
system at a rate proportional to the number of ions in the system. The sigmoid function is 
used to represent the growth and saturation of the streaming instability which causes 
bunching that increases the amplitude of the detected ion signal. The specific function 
that is normalized and plotted in figure 85 is 

Envelope(t) = ^ /0.0004^ Equation 4.5 

fl + exnf - ^ ~ 0.0001)/ Y| 

[ l+exp ^ /0.00006 J J 

Note the 100 ps time shift of the sigmoid to align the growth with the time axis of the 
data. Based on this model, the growth and saturation of the instability happens on a much 
faster time scale (60 ps) than the loss of ions from the system (400 ps). The growth of 
the instability in these data may be artificially accelerated, however, by the injection of 
the ion slug (an artificially large perturbation) at ion-source shut-off. 

Due to this different method of evaluation of confinement time in the presence of the 
instability at these lower pressures, these confinement time data are not included in the 
earlier presentations of confinement time data. There is also no good measurement of the 
two grid confinement time at this lowest pressure for purposes of comparison. If the 
trend at the lowest pressures is extrapolated from figure 59, the 2-grid confinement would 
be approximately 135ps at this pressure - a factor of 3 lower than the measured 400ps 
multi-grid time. 
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The specific conditions and mechanism for this instability are investigated in depth in the 
Ph.D. thesis of Tom McGuire. Experimentally, in the grid geometry tested, the instability 
was seen with a number of different confining potential structures. The best signal (most 
bunching) at all pressures was achieved with grids at the following potentials: [-5000, - 
2400, -830, -14, 0]. Double bunching was never seen in experiment although McGuire 
sees it in computational modeling. The detailed investigation of this phenomenon is 
beyond the scope of this work. 
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5. Conclusions and Future Work 


Computational and experimental investigations into improving ion confinement in 
inertial electrostatic confinement systems were conducted. The baseline 2-grid system 
was compared to a multi-grid system and it was found that the shaping of the electrostatic 
background field could improve ion confinement in the high vacuum operating regime. 
The PIC predicted ion bunching mode which is believed to be a saturated two-stream 
instability was also seen for the first time in an IEC experiment. Computational 
predictions of 2-grid confinement using OOPIC were quite poor, but modeling of the 
pressure sensitivity of multi-grid confinement was matched quite well in experiment 
down to pressures at which the two-stream instability was seen to develop. At these 
pressures (below ~4e-5 mbar), the 1/e confinement was difficult to clearly distill from the 
data using the same techniques due to the presence of the instability. An envelope 
amplitude technique was used to show that a multi-grid Argon ion confinement test at 
1.9e-6 mbar had a 1/e confinement time of approximately 400ps, or approximately 33 
passes compared to the two grid baseline in which the 1/e confinement time was shown 
to be approximately 10 passes at best. 

5. 1. Contributions of this thesis 

This thesis reports on the modeling, design, construction, and testing of the first 
published multi-grid IEC ion confinement experiment. The general scaling for non- 
neutral IEC systems was derived. A useful non-dimensional number for non-neutral IEC 
systems was proposed. The design and construction of the first multi-grid IEC 
experiment was reported. A 3-dimensional, multi-grid semi-analytic potential model was 
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described. The first direct measurements of ion confinement times in a multi-grid IEC 
device were reported. Computational experiments of particle-in-cell models created in 
OOPIC were compared to measurements of a physical experiment, and the potential for 
multi-grid IEC devices to improve ion confinement times was shown both 
computationally and in experiment. Detection of the computationally predicted two- 
stream instability at the ion bounce frequency was also reported. 

5.2. Suggestions 

The current experimental hardware could be used to further these investigations into the 
confinement properties of multi-grid IEC systems at HV and UHV pressures more 
effectively if an improved ion source with a much higher pressure in the ionization region 
than in the rest of the chamber is developed. This type of source could decouple the ion 
signal strength from the background density which could allow for much more accurate 
measurements at much lower pressures. It is also desirable to redesign the ion source to 
eliminate the slug of ions that is injected when the source is shut down. 

It would be interesting to experimentally explore the conditions for the onset of the 
instability. This type of investigation would be greatly informed by continued theoretical 
research such as is embodied in the work of McGuire [357]. 

Because of the scaling of this type of non-neutral IEC device, improved power outputs 
could be achieved by shrinking the size of the device. The limits to this shrinkage would 
be imposed by surface flashover considerations, or possibly by cold Fowler-Norheim 
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emission of electrons from the cathode resulting in rapid loss of efficiency at very small 
scale. Such analysis is left to future students. 
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Appendix: Code 

This appendix contains some of critical pieces of code that were generated in order to 
produce the computational results reported in this thesis. Code is written for MATLAB 
Version 14 and OOPIC Pro Version 1.0. For compactness of this document, only one 
version of the OOPIC model is included (GridDiagnostics5.inp). 
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Matlab Code 


For reference I have only included the MATLAB code which was directly used in the 
Simulated Annealing design routine. Elements of this code are used in other scripts to 
plot parameters of interest, but the included package is self-sufficient for the SA design 
process. 


Annealit2thengrad.m 

% Simulated Annealing design script for IEC fusion experiment 
% Carl Dietrich 5/31/05 

% It is assumed that there are 8 longitude lines and 4 latitude lines of 
% even spacing (45 and 36 degree separations), and that the 5th, outermost grid is 
% grounded at a radius of 0.25 m. 

% Configuration = [Ro Rinjector Rgridl Rgrid2 Rgrid3 Rgrid4 Potgridl Potgrid2 
% Potgrid3 Potgrid4] ; 

clear all 
close all 

Ro = .01; 

Rinjector = .18; 

Rgridl = .05; 

Rgrid2 = .07; 

Rgrid3 = . 1 ; 

Rgrid4 = .14; 

Potgridl = -7000; 

Potgrid2 = -10000; 

Potgrid3 = -4000; 

Potgrid4 = -1000; 

Configuration = [Ro Rinjector Rgridl Rgrid2 Rgrid3 Rgrid4 Potgridl Potgrid2 Potgrid3 
Potgrid4] ; 

%[xbest,Ebest,xhist]=SA(xo,file_eval,file_perturb, options); 

% options algorithm option flags. Uses defaults, [ ], if left blank 
% (1) To - initial system temperature - automatically determined if 

% left blank ([]). To should be set such that the expression 

% exp(-E(xo)/To)>0.99 is true, i.e. the initial system is "melted" 

% (2) Cooling Schedule: linear=l, exponential [2] 

% (3) dT Temp, increment, e.g. [dT=0.9] for exp. cooling Tk=dT A k*To, 

% abs. temperature increment for linear cooling (Tk+l=Tk-dT); 

% (4) neq = equilibrium condition, e.g. number of rearrangements 
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% attempted to reach equilibrium at a given temperature, neq=[5] 

% (5) frozen condition - sets up SA exit criterion 

% nfrozen = non-integer, e.g. 0. 1 SA interprets this numbers as Tmin, 

% the minimum temperature below which the system is frozen. 

% nfrozen = integer ,e.g. 1,2.. SA interprets this as # of successive 

% temperatures for which the number of desired acceptances defined 

% under options(4) is not achieved, default: nfrozen=[3] 

% (6) set to 1 to display diagnostic messages (=[1]) 

% (7) set to 1 to plot progress during annealing (=[0]) 

To=le2; options(l)=To; 
schedule=2; options(2)=schedule; 
dT=.l; options(3)=dT; 
neq=100; options(4)=neq; 
nfrozen=0.001 ; options(5)=nfrozen; 
diagnostics=l ; options(6)=diagnostics; 
plotflag=l; options(7)=plotflag; 

rand('state',sum( 100*clock)); 
tic 

[BestConfig,Ebest,xhist]=SAmodCarl(Configuration,'griddesigneval2','PerturbConfig2',o 

ptions) 

toe 

[a b]=size(BestConfig); 

PlotConfig2(BestConfig(a,:)); 
title('Modified SA Results'); 

coreacceptance = .5; 
mingap = .02; 

LB(1)=0; 

UB( 1 )=BestConfig(a,3)*sin(pi/8)*coreacceptance; %max Ro is the first grid radius times 
the acceptance 

LB(2)=BestConfig(a,5); % lower bound of Rinj is 3rd grid radius 
UB(2)=.25; 

LB(3)=BestConfig(a,3); %1.5*mingap; 

UB(3)=(BestConfig(a,4)+BestConfig(a,3))/2 - mingap/2; 

LB(4)=UB(3)+mingap; 

UB(4)=(BestConfig(a,5)+BestConfig(a,4))/2 - mingap/2; 

LB(5)=UB(4)+mingap; 

UB(5)=(BestConfig(a,6)+BestConfig(a,5))/2 - mingap/2; 

LB(6)=UB(3)+mingap; 

UB(6)=.25-mingap; 

LB(7)=- 10000; 

UB(7)=0; 

LB(8)=- 10000; 

UB(8)=0; 

LB(9)=- 10000; 
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UB(9)=0; 

LB( 10)=- 10000; 

UB(10)=0; 

tic 

[NewBestConfig,LowestEnergy] = 

fmincon('griddesigneval2',BestConfig(a, :) , [],[],[],[] ,LB ,UB ) 
toe 

PlotConfig2(NewBestConfig) ; 
title('FMINCON Gradient Search Results'); 

% Save results to file: 

dumpfilename = strcat('S A_FMINCON_',datestr(now,'mmm-dd-HHMM')) ; 
fid=fopen(dumpfilename, ' w') ; 
if fid==-l 

error = 'Unable to open dump file!' 
end 

fid=fopen('filename' , ' w') 

fprintf(fid,'%s \n',num2str(NewBestConfig)); 

fclose(fid); 
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Griddesigneval2.m 


% Energy evaluation function for simulated annealing design of multi-grid 
% IEC fusion reactor experiment 
% Carl Dietrich 5/29/05 

% Input is an array 

% Configuration = [Ro Rinjector Rgridl Rgrid2 Rgrid3 Rgrid4 Potgridl Potgrid2 
% Potgrid3 Potgrid4] ; 

% It is assumed that there are 8 longitude lines and 4 latitude lines of 
% even spacing (45 degree separations), and that the 5th, outermost grid is 
% grounded at a radius of 0.2 m. 

function Energy = griddesigneval(Configuration) 

global rx pot_beam pot_prime_beam pot_dubprime_beam K inj_pot 

e=1.602e-19; 

epso=8.854e-12; 

c=3e8; 

% Input the desired grid structure: 

r = [Configuration^) Configuration^) Configuration^) Configuration(6) .25]; % grid 
radius vector 

if r(l)>r(2) I r(2)>r(3) I r(3)>r(4) I r(4)>r(5) 
disp('Waming: Invalid grid configuration!') 
r = r 

Energy = 10000 
return 
end 

Voltagedesired = [Configuration^) Configuration^) Configuration(9) Configuration 10) 

0 ]; 

n=max(size( Y oltagedesired)) ; 
solidityfudgefactor = 2.75e-5; 
for i=l:n-l 

deltaV(n-i) = -(Voltagedesired(n+l-i)-Voltagedesired(n-i)); 

Qins(n-i) = solidityfudgefactor*4*pi*epso*deltaV (n-i)/( l/r(n-i)- l/r(n+ 1 -i)) ; 
end 

Qins(n)=Qins(n- 1 ) ; 

Qdesired( 1 )=Qins( 1 ) ; 
for i=2:n 

Qdesired(i) = Qins(i)-Qins(i-1); 
end 

q = Qdesired./r. A 2;%[10.5e-8 -2.1e-7 7.5e-8 0]; % grid charge vector 

t = 4; % number of latitude wires 

1 = 8; % number of longitude wires 

thickness =.001;% wire thickness 

kmax = 32; 

mult = kmax/(2*l); 
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% Calculate coefficients: 

%tic 

[Abar, Bbar]=solve_coeff(r,q,t,l, thickness, mult); 

%toc 

%tic 

phi = pi/2; 
theta = 0; %pi/l; 
divs = 300; 
dr = max(r)/divs; 
rx=[0:dr:max(r)]; 

for i=l:(divs+l) 
if rx(i)<.01 

pot(i) = spherical_potential_norm_posm_bar(r, Abar,Bbar,.0 1 , theta, phi); 
else 

pot(i) = spherical_potential_norm_posm_bar(r, Abar, Bbar,rx(i), theta, phi) ; 
end 
end 

edge_ground = -pot(divs+l); 
pot = pot+edge_ground; 

for i=l:(divs+l) 
if i==l 

pot_prime_wire(i)=0; 
elseif i==divs+l 

pot_prime_wire(i)=(pot(i)-pot(i- 1 ))/dr ; 
else 

pot_prime_wire(i)=((pot(i) -pot(i- 1 ))/ dr+(pot(i+ 1 ) -pot(i))/ dr)/2 ; 
end 
end 

theta = pi/1; 
for i=l:(divs+l) 
if rx(i)<.01 

pot_beam(i) = spherical_potential_norm_posm_bar(r,Abar,Bbar,.0 1 , theta, phi); 
else 

pot_beam(i) = spherical_potential_norm_posm_bar(r, Abar,Bbar,rx(i) , theta, phi) ; 
end 
end 

pot_beam = pot_beam+edge_ground; 
for i=l:(divs+l) 
if i==l 

pot_prime_beam(i)=0 ; 
elseif i==divs+l 

pot_prime_beam(i)=(pot_beam(i)-pot_beam(i-l))/dr; 

else 
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pot_prime_beam(i)=((pot_beam(i)-pot_beam(i-l))/dr+(pot_beam(i+l)- 

pot_beam(i))/dr)/2; 

end 

end 

for i=l:(divs+l) 
if i==l 

pot_dubprime_beam(i) = 0; 
elseif i==divs+l 

pot_dubprime_beam(i)=(pot_prime_beam(i)-pot_prime_beam(i-l))/dr; 

else 

pot_dubprime_beam(i)=((pot_prime_beam(i)-pot_prime_beam(i- 
1 ))/dr+(pot_prime_beam(i+ 1 ) -pot_prime_beam(i) )/dr)/2 ; 

end 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

K=zeros(size(rx)); 

% Initial beam envelope conditions 
Ro=Configuration( 1 ) ; 
r_inj = Configuration(2); 
mi = 40*1.67e-27; % Argon 

inj_pot = spherical_potential_norm_posm_bar(r, Abar,Bbar,r_inj , pi/1, pi/2)+edge_ground; 
inj_pot_dtheta = 

spherical_potential_norm_posm_bar(r,Abar,Bbar,r_inj,2*pi/l,pi/2)+edge_ground; 

dphidtheta=(inj_pot_dtheta-inj_pot)/(pi/l); 

if dphidtheta<0 

disp('Reflection region is defocussing => bad design') 

Energy = 10000 
return 
end 

vtherm=340; 

Rprimeo=0; 

Rfinal = 0; 

Rprimefinal = 0; 
no=2el3; 

YO = [Ro Rprimeo]; 

OPTIONS = ODESET ('RelTol' , 1 e-4, ' AbsTol' , 1 e- 8) ; 
for i=l:max(size(rx)) 

vi(i)=sqrt(2*e*(abs(pot_beam(i)-inj_pot))/mi)+vtherm; 
gamma(i)= l/sqrt(l-vi(i) A 2/c A 2); 

K(i)=e A 2*no*Ro A 2/(2*epso*mi*vi(i) A 2*gamma(i) A 3); % perveance 
end 

if max(isnan(K)) 
disp('K is NaN!') 

Energy = 10000 
return 
end 
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iter = 0; 

maxiters = 1000; 
acceptance = .8; 

while itercmaxiters % & beam envelope is inside acceptance window 
iter = iter+1; 
no = l.l*no; 

K = K.*1.2; % since perveance is proportional to density... 

[X Y]=ode45('ode45beamfunc',[0 (r_inj-dr)],YO, OPTIONS); 

[n m]=size(Y); 

[maxR, imaxR]=max(abs(Y(:,l))); 

%Checks: 

if maxR > X(imaxR)*sin(pi/8)* acceptance & X(imaxR)>=r(l) 
break 
end 

for ib=l:4 

if r_inj>r(ib) & interpl(X, Y(:,l), r(ib)*cos(pi/8)) > acceptance *r(ib)*sin(pi/8) 
disp('good break') 
break 
end 
end 
end 


no = no; 
Ro = Ro; 


normalizing_factor = (lel4) A 2*(.01) A 3; % Approximate density and core size from 2D 
OOPIC sims 

Energy = normalizing_factor/(no A 2*Ro A 3) %*V_c A 4) 
if max(max(isnan(Y))) 

disp('Waming: Solution NaN!') 

Energy = 10000 
return 
end 

if min(Y(:,l))<0 

disp('Waming: beam envelope crosses R=0!') 

Energy = 10000 
end 

if max(imag( Y ( : , 1 ))) 

disp('Beam envelope has imaginary component!') 

Energy = 10000 
end 

if iter>=maxiters 

disp('Maximum number of iterations has been exceeded!') 

Energy = 10000 
end 
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%toc 
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PerturbConfig2.m 


% Configuration perterbation function for SA design of a multi-grid 
% IEC fusion reactor experiment 
% Carl Dietrich 5/31/05 

% Input is an array 

% Configuration = [Ro Rinjector Rgridl Rgrid2 Rgrid3 Rgrid4 Potgridl Potgrid2 
% Potgrid3 Potgrid4] ; 

% It is assumed that there are 8 longitude lines and 4 latitude lines of 
% even spacing (45 degree separations), and that the 5th, outermost grid is 
% grounded at a radius of 0.25 m. 

function [NewConfig] = PerturbConfig2(Configuration) 

% Start with old configuration: 

NewConfig = Configuration; 

% Pick the two degrees of freedom to be perturbed 
firstchoice = ceil(10*rand); 
secondchoice = ceil(10*rand); 
while firstchoice==l 

firstchoice = ceil(10*rand); 
end 

while firstchoice == secondchoice I secondchoice==l % ensure that there are 2 perturbed 
DOFs 

secondchoice = ceil(10*rand); 
end 

% If it is necessary to move the grids, move them first (inside out)... 
mingap = .02; % minimum radial gap between grids... 

if firstchoice == 3 I secondchoice == 3% if we are perturbing the 1st grid 
% if Configuration^) -mingap <.25/3 

% maxsize=Configuration(4) -mingap; % can't move beyond the 2nd grid 

% if maxsize<minsize 

% maxsize=Configuration(3); 

% end 

% else 

% maxsize=.25/3; % or l/3rd the distance (arbitrary) 

% end 

minsize=.03; % minimum grid radius is 3 cm 

maxsize=Configuration(4)-mingap; % can't move beyond the 2nd grid 
NewConfig(3)= minsize+rand*(maxsize-minsize); 
if NewConfig(3)< minsize 
NewConfig(3)=minsize; 
elseif NewConfig(3)>maxsize 
N ewConfig(3 )=maxsize ; 
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end 

R_Grid_l = NewConfig(3) 
end 

if firstchoice == 4 I secondchoice == 4% if we are perturbing the 2nd grid 
% if Configuration(5)-mingap<.25/2 

% maxsize=Configuration(5)-mingap; % can't move beyond the 3rd grid 
% else 

% maxsize=.25/2; % or 1/2 the distance (arbitrary) 

% end 

maxsize=Configuration(5)-mingap; % can't move beyond the 3rd grid 
minsize=NewConfig(3)+mingap; % minimum grid radius is grid 1 + the min gap 
NewConfig(4)= minsize+rand*(maxsize-minsize); 

R_Grid_2 = NewConfig(4) 
end 

if firstchoice == 5 I secondchoice == 5% if we are perturbing the 3rd grid 
% if Configuration(6)-mingap<.25*.75 

% maxsize=Configuration(6)-mingap; % can't move beyond the 4th grid 
% else 

% maxsize=.25*.75; % or 3/4 the distance (arbitrary) 

% end 

maxsize=Configuration(6)-mingap; % can't move beyond the 4th grid 
minsize=NewConfig(4)+mingap; % minimum grid radius is grid2 +mingap 
NewConfig(5)= minsize+rand*(maxsize-minsize); 

R_Grid_3 = NewConfig(5) 
end 

if firstchoice == 6 I secondchoice == 6% if we are perturbing the 4th grid 
maxsize=.25-mingap; % max size is mingap from the anode 

minsize=NewConfig(5)+mingap; % minimum grid radius is grid3 plus mingap 
NewConfig(6)= minsize+rand*(maxsize-minsize); 

R_Grid_4 = NewConfig(6) 
end 

% Now pick focal point and injection point if necessary: 

% if firstchoice == 1 I secondchoice == 1% if we are perturbing the focal size 
% % maximum size of Ro (focus) is assumed to be half of the sine of the 

% % pi/8 times the inner grid radius, Rgridl 

% maxsize=0.5 *NewConfig(3) *sin(pi/8) ; 

% minsize=.003; % implicit assumption that the focal diameter > 6mm 
% NewConfig(l)= minsize+rand*(maxsize-minsize); 

% Ro = NewConfig(l) 

% end 

if firstchoice == 2 I secondchoice == 2% if perturbing the injector location 
% Assume the injector must be located outside the third grid and 
% inside the outermost anode grid (which has Ranode=0.25) 
maxsize = 0.25; 
min size = NewConfig(5); 

NewConfig(2)= minsize+rand*(maxsize-minsize); 

R_injector = NewConfig(2) 
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end 


% Now perturb the potentials if necessary: 

if firstchoice == 7 I secondchoice == 7% if we are perturbing the 1st pot 
maxpot=-5000; 
minpot=- 10000; 

NewConfig(7)= minpot+rand*(maxpot-minpot) ; 

Pot_Grid_l = NewConfig(7) 
end 

if firstchoice == 8 I secondchoice == 8% if we are perturbing the 1st pot 
maxpot=-5000; 
minpot=- 10000; 

NewConfig(8)= minpot+rand*(maxpot-minpot) ; 

Pot_Grid_2 = NewConfig(8) 
end 

if firstchoice == 9 I secondchoice == 9% if we are perturbing the 1st pot 
maxpot=0; 
minpot=-5000; 

NewConfig(9)= minpot+rand*(maxpot-minpot) ; 

Pot_Grid_3 = NewConfig(9) 
end 

if firstchoice ==101 secondchoice == 10% if we are perturbing the 1st pot 
maxpot=0; 
minpot=-5000; 

N ewConfig( 10)= minpot+rand* (maxpot-minpot) ; 

Pot_Grid_4 = NewConfig(lO) 
end 


■ Never perturb focal point (always constant fraction of inner grid radius) 
% maximum size of Ro (focus) is assumed to be half of the sine of the 
% pi/8 times the inner grid radius, Rgridl 
maxsize=0.5*NewConfig(3)*sin(pi/8); 

%minsize=.003; % implicit assumption that the focal diameter > 6mm 
NewConfig(l)= maxsize; %minsize4*and*(maxsize-minsize); 

Ro = NewConfig(l); 
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PlotConfig2.m 


% Energy evaluation function for simulated annealing design of multi-grid 
% IEC fusion reactor experiment 
% Carl Dietrich 5/29/05 

% Input is an array 

% Configuration = [Ro Rinjector Rgridl Rgrid2 Rgrid3 Rgrid4 Potgridl Potgrid2 
% Potgrid3 Potgrid4] ; 

% It is assumed that there are 8 longitude lines and 4 latitude lines of 
% even spacing (45 degree separations), and that the 5th, outermost grid is 
% grounded at a radius of 0.2 m. 

function Energy = PlotConfig2(Configuration) 

global rx pot_beam pot_prime_beam pot_dubprime_beam K inj_pot 

e=1.602e-19; 

epso=8.854e-12; 

c=3e8; 

% Input the desired grid structure: 

r = [Configuration^) Configuration^) Configuration^) Configuration(6) .25]; % grid 
radius vector 

if r(l)>r(2) I r(2)>r(3) I r(3)>r(4) I r(4)>r(5) 
disp('Invalid grid configuration!') 

Energy = 10000 
return 
end 

Voltagedesired = [Configuration^) Configuration^) Configuration(9) Configuration(lO) 

0 ]; 

n=max(size( Y oltagedesired)) ; 
solidityfudgefactor = 2.75e-5; 
for i=l:n-l 

deltaV(n-i) = -(Voltagedesired(n+l-i)-Voltagedesired(n-i)); 

Qins(n-i) = solidityfudgefactor*4*pi*epso*deltaV(n-i)/( 1 /r(n-i)- 1 /r(n+ 1 -i)); 
end 

Qins(n)=Qins(n- 1 ) ; 

Qdesired( 1 )=Qins( 1 ) ; 
for i=2:n 

Qdesired(i) = Qins(i)-Qins(i-1); 
end 

q = Qdesired./r. A 2;%[10.5e-8 -2.1e-7 7.5e-8 0]; % grid charge vector 

t = 4; % number of latitude wires 

1 = 8; % number of longitude wires 

thickness =.001;% wire thickness 

kmax = 32; 

mult = kmax/(2*l); 
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% Calculate coefficients: 

%tic 

[Abar, Bbar]=solve_coeff(r,q,t,l, thickness, mult); 

%toc 

%tic 

phi = pi/2; 
theta = 0; %pi/l; 
divs = 300; 
dr = max(r)/divs; 
rx=[0:dr:max(r)]; 

for i=l:(divs+l) 
if rx(i)<.01 

pot(i) = spherical_potential_norm_posm_bar(r, Abar,Bbar,.0 1 , theta, phi); 
else 

pot(i) = spherical_potential_norm_posm_bar(r, Abar, Bbar,rx(i), theta, phi) ; 
end 
end 

edge_ground = -pot(divs+l); 
pot = pot+edge_ground; 

for i=l:(divs+l) 
if i==l 

pot_prime_wire(i)=0; 
elseif i==divs+l 

pot_prime_wire(i)=(pot(i)-pot(i- 1 ))/dr ; 
else 

pot_prime_wire(i)=((pot(i) -pot(i- 1 ))/dr+(pot(i+ 1 ) -pot(i))/dr)/2 ; 
end 
end 
figure; 

subplot(2,l,l); 
plot(rx,pot,'b'); 
hold on; 
theta = pi/1; 
for i=l:(divs+l) 
if rx(i)<.01 

pot_beam(i) = spherical_potential_norm_posm_bar(r, Abar,Bbar, .0 1 , theta, phi) ; 
else 

pot_beam(i) = spherical_potential_norm_posm_bar(r, Abar, Bbar,rx(i), theta, phi); 
end 
end 

pot_beam = pot_beam+edge_ground; 
for i=l:(divs+l) 
if i==l 

pot_prime_beam(i)=0 ; 
elseif i==divs+l 

pot_prime_beam(i)=(pot_beam(i)-pot_beam(i-l))/dr; 
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else 

pot_prime_beam(i)=((pot_beam(i)-pot_beam(i-l))/dr+(pot_beam(i+l)- 

pot_beam(i))/dr)/2; 

end 

end 

for i=l:(divs+l) 
if i==l 

pot_dubprime_beam(i) = 0; 
elseif i==divs+l 

pot_dubprime_beam(i)=(pot_prime_beam(i)-pot_prime_beam(i-l))/dr; 

else 

pot_dubprime_beam(i)=((pot_prime_beam(i)-pot_prime_beam(i- 
1 ))/dr+(pot_prime_beam(i+ 1 ) -pot_prime_beam(i))/dr)/2 ; 

end 

end 

plot(rx,pot_beam,'r'); 
title('Potential v. radius'); 

legend('Theta = 0 (on grid)', 'Theta = pi/1 (on beampath)'); 

xlabel('r(m)'); 

ylabel('Potential (V)'); 

drawnow; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

K=zeros(size(rx)) ; 

% Initial beam envelope conditions 
Ro=Configuration( 1 ) ; 
r_inj = Configuration^); 
mi = 40* 1 .67e-27 ; % Argon 

inj_pot = spherical_potential_norm_posm_bar(r, Abar,Bbar,r_inj , pi/1, pi/2)+edge_ground; 
inj_pot_dtheta = 

spherical_potential_norm_posm_bar(r, Abar,Bbar,r_inj ,pi/l+pi/ 1 000,pi/2)+edge_ground; 
dphidtheta=(inj_pot_dtheta-inj_pot)/(pi/ 1 000) ; 
if dphidthetacO 

disp('Reflection region is defocussing => bad design') 

Energy = 10000 
return 
end 

vtherm=340; 

Rprimeo=0; 

Rprimefinal=0; 

Rfinal = 0; 
no=2el3; 

YO = [Ro Rprimeo]; 

OPTIONS = ODESET ('RelTol' , 1 e-4, ' AbsTol' , 1 e- 8) ; 
for i=l:max(size(rx)) 

vi(i)=sqrt(2*e*(abs(pot_beam(i)-inj_pot))/mi)+vtherm; 
gamma(i)= l/sqrt(l-vi(i) A 2/c A 2); 

K(i)=e A 2*no*Ro A 2/(2*epso*mi*vi(i) A 2*gamma(i) A 3); % perveance 
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end 

if max(isnan(K)) 
disp('K is NaN!') 

Energy = 10000 
return 
end 

% X=0; 

% Xback=0; 

% maxR=0; 

% imaxR=l; 

% maxRbaek=0; 

% imaxRback=l; 

% acceptance = .8; 

% while Rfinal<Ro & maxR<=acceptance*X(imaxR)*tan(pi/8)& 
maxRback<=acceptance*Xback(imaxRback)*tan(pi/8) 

% no = 2*no; 

% K = K.*2; % since perveance is proportional to density... 

% [X Y] =ode45 ('ode45beamfunc ' , [0 (rjnj -dr) ],YO, OPTIONS); 

% [n m]=size(Y); 

% Roback = Y(n,l); 

% Rprimeback = Y(n,2); 

% [dphidx dphidy]=gradfield(r_inj-dr,pi/l+atan(Roback/(r_inj- 
dr)) , Abar,Bbar,r,edge_ground) ; 

% Omega = atan(dphidy/dphidx); 

% omega_one = atan(Rprimeback); 

% omega_two = 2*Omega - omega_one; 

% Rprimeretum = tan(omega_two); % mirror reflection off of equipotential 
% [Xback Yback]=ode45('ode45beamfunc',[(r_inj-dr) 0], [Roback 
Rprimeretum] , OPTIONS) ; 

% [n m]=size(Yback); 

% Rfinal = Yback(n, 1 ) ; 

% [maxR, imaxR]=max(Y(:,l)); 

% [maxRback, imaxRback]=max(Yback(:, 1 )); 

% end 
iter = 0; 
maxiters = 30; 
acceptance = .8; 

while iter<maxiters % & beam envelope is inside acceptance window 
iter = iter+1; 
no = l.l*no; 

K = K.*l .2; % since perveance is proportional to density... 

[X Y]=ode45('ode45beamfunc',[0 (r_inj-dr)],YO, OPTIONS); 

[n m]=size(Y); 

[maxR, imaxR]=max(abs(Y(:,l))); 

%Checks: 

if (maxR > X(imaxR)*sin(pi/8)*acceptance & X(imaxR)>=r(l)) 
break 
end 
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for ib=l:3 

if interpl(X, Y(:,l), r(ib)*cos(pi/8)) > acceptance *r(ib)*sin(pi/8) 
break 
end 
end 
end 
%toc 
no = no 
Ro = Ro 

%Rprimefinal = Rprimefinal 

normalizing_factor = (lel4) A 2*(.01) A 3; % Approximate density and core size from 2D 
OOPIC sims 

Energy = normalizing_factor/(no A 2*Ro A 3) %*V_c A 4) 
if max(max(isnan(Y))) 

disp('Waming: Solution NaN!') 

Energy = 10000 
return 
end 

if min(Y(:,l))<0 

disp('Waming: beam envelope crosses R=0!') 

Energy = 10000 
end 

if max(imag(Y(:,l))) 

disp('Beam envelope has imaginary component!') 

Energy = 10000 
end 

if iter>=maxiters 

disp('Maximum number of iterations has been exceeded!') 

Energy = 10000 
end 


subplot(2,l,2); 
plot(r_inj,0,'g*'); 
hold on; 

plot(X,Y(:,l),'g'); 

%plot(Xback, Yback( : , 1 ),'b') ; 
for i=l:max(size(r)) 

plot(r(i) *cos(pi/8) ,r(i) * sin(pi/8) , 'r * ') ; 
end 

plot(rx. *cos(pi/8) ,rx. *sin(pi/8) ,'r') ; 
labelone=sprintf('Gridline') ; 

labeltwo = sprintf('Outward Beam Envelope for n_c=%3.2E m-3, r_o=%3.2E m',no,Ro); 
%labelthree = sprintf('Inward Beam Envelope for n_o=%3.2E m-3, r_o=%3.2E 
m',no*vi( 1 )/vtherm*(Ro/Roback) A 2,Roback); 
legend('Injection Point', labeltwo, labelone) ; 
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SAmodCarl.m 


% Note: this code is based upon and nearly identical to SA.m developed by Prof. Olivier 
de Week 

function [xbest,Ebest,xhist]=SAmodCarl(xo,file_eval,file_perturb, options); 

% [xbest,Ebest,xhist]=SA(xo,file_eval,file_perturb, options); 

% Single Objective Simulated Annealing (SA) Algorithm 

% This function is a generic implementation of the single objective 
% Simulated Annealing (SA) algorithm first proposed by Kirkpatrick, 

% Gelatt and Vecchi. The algorithm tries to improve upon an initial 
% configuration, xo, by evaluating perturbed configurations. When the 
% system reaches the "frozen" state, the algorithm stops and the best 
% configuration and search history are returned. The user can choose 
% from one of two cooling schedules: linear or exponential. 

% Input: 

% xo initial configuration of the system (a row vector) 

% file_eval file name (character string) of configuration evaluator; 

% assumes that E='file_eval'(x) is a legitimate function 

% call; set up function such that (scalar) output E will be 

% minimized. 

% file_perturb file name (character string) of configuration perturbator; 

% assumes that xp='fname_perturb'(x) is a legitimate function 

% call. This function creates a "neighboring" configuration. 

% options algorithm option flags. Uses defaults, [ ], if left blank 
% (1) To - initial system temperature - automatically determined if 

% left blank ([]). To should be set such that the expression 

% exp(-E(xo)/To)>0.99 is true, i.e. the initial system is "melted" 

% (2) Cooling Schedule: linear=l, exponential [2] 

% (3) dT Temp, increment, e.g. [dT=0.9] for exp. cooling Tk=dT A k*To, 

% abs. temperature increment for linear cooling (Tk+l=Tk-dT); 

% (4) neq = equilibrium condition, e.g. number of rearrangements 

% attempted to reach equilibrium at a given temperature, neq=[5] 

% (5) frozen condition - sets up SA exit criterion 

% nfrozen = non-integer, e.g. 0. 1 SA interprets this numbers as Tmin, 

% the minimum temperature below which the system is frozen. 

% nfrozen = integer ,e.g. 1,2.. SA interprets this as # of successive 

% temperatures for which the number of desired acceptances defined 

% under options(4) is not achieved, default: nffozen=[3] 

% (6) set to 1 to display diagnostic messages (=[1]) 

% (7) set to 1 to plot progress during annealing (=[0]) 

% 

% Output: 

% xbest Best configuration(s) found during execution - row vector(s) 

% Ebest Energy of best configuration(s) (lowest energy state(s) found) 
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o xhist 
o .iter 


.E 

.T 

.k 

.C 

.S 


structure containing the convergence history 
Iteration number (number of times file_eval was called) 
current configuration at that iteration 
current system energy at that iteration 
current system temperature at that iteration 
temperature step index k 
specific heat at the k-th temperature 
entropy at the the k-th temperature 


% .Tnow temperature at the k-th temperature step 

% 

% User Manual (article): SA.pdf 

% 

% Demos: SAdemoO - four atom placement problem 

% SAdemol - demo of SA on MATLAB peaks function 

% SAdemo2 - demo of SA for Travelling Salesman Problem (TSP) 

% SAdemo3 - demo of SA for structural topology optimization 

% SAdemo4 - demo of SA for telescope array placement problem 


% dWo,(c) MIT 2004 


% Ref: Kirkpatrick, S., Gelatt Jr., C.D. and Yecchi, M.P., "Optimization 
% by Simulated Annealing", Science, Vol. 220, Number 4598, pp. 671-680, May 
% 1983 


% dump file added by Carl Dietrich 6/1/05 
dumpfilename = datestr(now,'mmm-dd-HHMM') 
fid=fopen(dumpfilename, ' w') ; 
if fid==-l 

error = 'Unable to open dump file!' 
end 

%%%%% 


%check input 
if -isempty(options) 

To=options(l); 
schedule=options(2) ; 
dT=options(3); 
neq=options(4); 
nfrozen=options(5) ; 
diagnostics=options(6) ; 
plotflag=options(7) ; 
else 

% set all options to default 
% To - set initial system temperature 
eval(['Eo=' file_eval '(xo);']); 

To=abs(-Eo/log(0.99)); % set initial temperature such that probablility of 

% accepting an inferior solution is initially equal to 0.99 

schedule=2; 
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dT=0.9; % this is the ratio dT=(T_i+l/T_i) for geometrical cooling 
neq=5; % number of rearrangements accepted at a given T 
nffozen=3; % if neq hasn't been reached at nffozen successive 
% temperatures the system is considered frozen and the SA exits 
diagnostics=l; % display messages 
plotflag=0; %plot convergence 
end 
% 

nmax=neq*round(sqrt(max(size(xo)))); % nmax - maximum number of steps at one 
temperature, while 

% trying to establish thermal equilibrium 

if nfrozen==round(nfrozen) 

% nfrozen is integer - look for nfrozen successive temperatures without 
% neq acceptances 
Tmin=0; 
else 

T min=nfrozen ; nfrozen=3 ; 
end 


% Step 1 - Show initial configuration 

if diagnostics==l 

disp('Initial configuration: ') 

xo 

end 

% Step 2 - Evaluate initial configuration 
eval(['Eo=' file_eval '(xo);']); 
counter=l; 

xnow=xo; Enow=Eo; nnow=l; 
xhist(nno w) . iter=counter ; 

xhist(nno w) . x=xo ; 
xhis t(nno w ) . E=Eno w ; 
xhist(nnow).T=To; 

% still need to add .S current entropy at that iteration 
xbest=xnow; 

fprintf(fid,' % s \n' ,num2str(xbest)) ; 

Ebest=Enow; 

Tnow=To; 
if diagnostics==l 

disp(['Energy of initial configuration Eo: ' num2str(Eo)]) 
end 

if plotflag 
figure(99) 
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semilogy(counter,Enow,'k*'); %plot(counter,Enow,'k*'); 
hold on 

semilogy(counter,Enow,'mo') %plot(counter,Enow,'mo') 
xlabel('Iteration Number') 
ylabel('System Energy') 

legend('current configuration', 'new best configuration') 
title('SA convergence history') 
lastbest=counter; 
drawnow 
end 

frozen=0; % exit flag for SA 

naccept=l; % number of accepted configurations since last temperature change 
Tlast= 1 ; % counter index of last temperature change 
k= 1 ; % first temperature step 

ET=[]; % vector of energies at constant system temperature 

% start annealing 

while (frozen<nfrozen)&(Tnow>Tmin) 

%Step 3 - Perturb xnow to obtain a neighboring solution 
if diagnostics 

disp(['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' Perturbing 
configuration']) 
end 

eval(['xp=' file_perturb '(xnow);']); 

%Step 4 - Evaluate perturbed solution 
eval(['Ep=' file_eval '(xp);']) 
counter=counter+ 1 ; 

%Step 5 - Metropolis Step 

dE=Ep-Enow; % difference in system energy 

PdE=exp(-dE/Tnow); 

if diagnostics 

disp(['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' P(dE)= ' num2str(PdE)]) 
end 

%Step 6 - Acceptance of new solution 

if dE<=0 % energy of perturbed solution is lower , automatically accept 
nnow=nnow+l; 
xnow=xp; Enow=Ep; 
xhist(nnow) . iter=counter ; 
xhist(nno w) . x=xp ; 
xhist(nnow).E=Ep; 
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xhist(nnow).T=Tnow; 
naccept=naccept+ 1 ; 
if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' Automatically accept 
better configuration (downhill)']) 
end 

else 

% energy of perturbed configuration is higher, but might still accept it 
randomnumberO 1 =rand; 
if PdE>randomnumber01 
nnow=nnow+l; 
xnow=xp; Enow=Ep; 
xhist(nno w) . iter=counter ; 
xhist(nnow) .x=xp ; 
xhist(nnow).E=Ep; 
xhist(nnow).T=Tnow; 
if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' Accepted inferior 
configuration (uphill)']) 
end 

else 

% keep current configuration 
xnow=xnow; 

Enow=Enow; 
if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' Kept the current 
configuration']) 
end 
end 
end 

ET=[ET; Enow]; 
if plotflag 
figure(99) 

semilogy(counter,Enow,'k*'); %plot(counter,Enow,'k*'); 

drawnow 

end 


if Enow<Ebest 

% found a new 'best' configuration 
Ebestlast=Ebest; 

Ebest=Enow; 

xbest=xnow; 

fprintf(fid,'%s \n',num2str(xbest)); 
if diagnostics 
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disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' This is a new best 
configuration']) 
end 

if plotflag 
figure(99) 

semilogy(counter,Enow,'mo') ; %plot(counter,Enow,'mo') ; 
semilogy([lastbest counter], [Ebestlast Enow],'m-'); %plot([lastbest 
counter], [Ebestlast Enow],'m-'); 
lastbest=counter; 
drawnow 
end 

elseif Enow==Ebest 
same=0; 

for ib=l :size(xbest, 1) 
if xbest(ib,:)==xnow 
if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' Found same best 
configuration']) 
end 

same=l ; 
end 
end 


if same ==0 
Ebestlast=Ebest; 

Ebest=Enow; 
xbest=[xbest ; xnow]; 
if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' Found another best 
configuration']) 
end 

if plotflag 
figure(99) 

semilogy(counter,Enow,'mo') ; %plot(counter,Enow,'mo') ; 
semilogy([lastbest counter], [Ebestlast Enow],'m-'); %plot([lastbest 
counter], [Ebestlast Enow],'m-'); 
lastbest=counter; 
drawnow 
end 
end 
end 

%Step 7 - Adjust system temperature 
Told=Tnow; 

if (naccept<neq)&(counter-Tlast)<nmax 
if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' Need to reach 
equilibrium at this temperature']) 
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end 

% continue at the same system temperature 
elseif (naccept<neq)&(counter-Tlast)>=nmax 
if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' System nearly frozen']) 
end 

Eavg=mean(ET); 

Evar=mean(ET. A 2) ; 

C=(Evar-Eavg A 2)/Tnow A 2; % specific heat 

S=log(nmax*length(unique(ET))/length(ET)); 

xhist(k).k=k; 

xhist(k).C=C; 

xhist(k).S=S; 

xhis t(k) .Tnow=Tnow; 


frozen=frozen+l; 

Tlast=counter; 

naccept=0; 


switch schedule 
case 1 

% linear cooling 
Tnow=Tnow-dT; 
if TnowcO 

frozen=nfrozen; %system temperature cannot go negative, exit 
end 
case 2 

% exponential cooling 
Tnow=dT*Tnow; 
case 3 

Tindex=Tindex+ 1 ; 
if Tindex>size(Tuser,l) 

frozen=nfrozen; % have run through entire user supplied cooling schedule 
else 

T now=T user(T index, 1 ) ; 

neq=Tuser(Tindex,2); 

end 

otherwise 

disp('Erroneous cooling schedule choice - option(2) - illegal') 


k=k+l; 
if plotflag 
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figure(98) 

hist(ET); Nh=hist(ET); Nh=max(Nh); 
hold on 

semilogx([Eavg Eavg]',[0 Nh+l],'k— ') %plot([Eavg Eavg]',[0 Nh+l],'k— ') 
text(Eavg, Nh+1, ['T=' num2str(Told,2)]) 
xlabel(Energy') 
ylabel('Occurences') 
drawnow 
end 
ET=[]; 


elseif (naccept==neq) 
if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' System reached 
equilibrium']) 
end 

Eavg=mean(ET); 

Evar=mean(ET. A 2) ; 

C=(Evar-Eavg A 2)/Tnow A 2; % specific heat 

S=log(nmax*length(unique(ET))/length(ET)); 

xhist(k).k=k; 

xhist(k).C=C; 

xhist(k).S=S; 

xhist(k) . Tnow=T now ; 


Tlast=counter; 

naccept=0; 

switch schedule 
case 1 

% linear cooling 
Tnow=Tnow-dT; 
if TnowcO 

frozen=nfrozen; %system temperature cannot go negative, exit 
end 
case 2 

% exponential cooling 
Tnow=dT*Tnow; 
case 3 

% user supplied cooling 
Tindex=Tindex+ 1 ; 
if Tindex>size(Tuser,l) 

frozen=nfrozen; %have run through entire user supplied cooling schedule 
else 

T now=T user(T index, 1 ) ; 
neq=Tuser(Tindex,2); 
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end 


otherwise 

disp('Erroneous cooling schedule choice - option(2) - illegal') 

end 


k=k+l; 

if plotflag 
figure(98) 

hist(ET); Nh=hist(ET); Nh=max(Nh); 
hold on 

semilogx([Eavg Eavg]',[0 Nh+l],'k~') %plot([Eavg Eavg]',[0 Nh+l],'k-') %([Eavg 
Eavg]',[0 Nh+l],'k— ') 

text(Eavg, Nh+1, ['T=' num2str(Told,2)]) 
xlabel('Energy') 
ylabel('Occurences') 
drawnow 
end 


ET=[]; 

end 

end %while (frozen<nfrozen)&(Tnow>tmin) 

fprintf(fid,' % s \n' ,num2str(xbest)) ; 
fclose(fid); 

% Reached end of SA 
if plotflag 
figure(97) 
k=k- 1 ; 
for ind=l:k 

S(ind)=xhist(ind).S; 

C(ind)=xhist(ind) . C ; 

T now(ind)=xhist(ind) .Tnow ; 
end 

plot([l:k],C,'bo') 
hold on 

plot([l:k],S,'ms') 

plot( [ 1 :k] ,log(Tnow) ,'kd') 

legend('C- specific heat' , 'S -entropy ' , 'ln(T) -temperature') 

xlabel('Temperature Step') 

title('Simulated Annealing Evolution') 

plot([l:k],C,'b-') 

plot([l:k],S,'m-') 

plot( [ 1 :k] ,log(Tnow) ,'k-') 
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drawnow 


end 

if diagnostics 

disp( ['Counter: ' num2str(counter) ' Temp: ' num2str(Tnow) ' System frozen, SA 
ended']) 

disp(['Best configuration: ']) 
xbest 

disp(['Lowest System Energy: ' num2str(Ebest) ]) 
end 
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Solve_coeff.m 


% Matlab function to calculate the coefficients of the potential 
% structure around multiple charged, concentric, spherical grids. 

% Written by: Carl Dietrich, 8/6/04 

% Inputs: r (radius vector of grids), q (charge/length on grids — same size 
% as r), t (number of latitude lines), 1 (number of longitude lines), 

% acc (the accuracy needed — typically .05). NOTE: the r vector must have 
% increasing radii with index and the charges in q correspond to the 
% respective r indices. 

% thickness is the wire thickness 
% mult = kmax/(2*l) 

% Outputs: A and B are 3d matrices of the coefficients for phi. 

% The first dimension specifies the region of validity in radius. 

% There will be n+1 regions. The second dimension is k and the 
% third dimension is m (2k+l). In the exact expression, k = inf 
% but in this model k is truncated to a finite integer determined 
% by the the required accuracy specified with the input 'acc'. 

function [Abar, Bbar]=solve_coeff(r,q,t,l, thickness, mult) 

kmax = 2*mult*l; 

stddevs_per_wirethickness = 2; 

numstddevs = 4; % 4 sigma = 0.9999366 of total 

% Check inputs: 
if size(r)~=size(q) 

error('Radius vector and charge vector must be the same length!'); 
end 
if t<l 

error('There must be at least one latitude line!'); 
end 
if 1<2 

error('There must be at least two longitude lines!') 
end 

if min(r)<thickness 

error('Geometry error: wires too big!'); 
end 

% Constants: 
epso = 8.854e-12; 

% Add in Ray's precalced integral: 
load IPkm 

% Find number of grids 
n = max(size(r)); 

% Calc effective wire thickness 

eff_thickness = numstddevs/stddevs_per_wirethickness*thickness; %pi/l*r(l) 
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% Calculate number of standard deviations per numerical integration 

%numstddevs = eff_thickness/thickness*stddevs_per_wirethickness 

divs=8; %50 %50 %100 % divisions per guassian (divs per eff_wire_thickness) 

% Initialize A and B 

Abar = zeros(n+l,kmax+l,kmax+l); 

Bbar = zeros(n+l,kmax+l,kmax+l); 

RHS = zeros(n+l,kmax+l,kmax+l); 
int = zeros(kmax+l,kmax+l); 

%wait = waitbar(0, 'Pre-calculating integrals for spherical potential solution...'); 

% Pre-Calc integral over phi for all k and m 
for k=0:kmax 

% waitbar(k/kmax,wait) ; 
for m=0:k 

int(k+l,m+l) = IPkm(k+l,m+l); 

end 

end 

%close(wait); 

%wait = waitbar(0, 'Calculating coefficients for spherical potential solution...'); 
for b=l:n 
for k=0:kmax 

%waitbar(k/kmax*b/n,wait) ; 

%if k==0 I rem(k,2)==0 % assumes north-south symmetry 
for m = 0:k 
RHS(b,k+l,m+l) = 0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% loop over longitude lines: 

sumint=0; 

for j= 1:1 

% sum of numerical integrals over cosmtheta 
inttwo=0; 

theta=2*pi*j/l-eff_thickness/r(b)/2; 
dtheta = eff_thickness/r(b)/divs; 
for i=l:divs 

inttwo=inttwo+gsinc(theta- 

2*pi*j/l,eff_thickness,r(b),numstddevs)*cos(m*theta)*dtheta; 
theta=theta+dtheta ; 
end 

sumint=sumint+inttwo ; 
end 

RHS(b,k+l,m+l) = RHS(b,k+l,m+l)+int(k+l,m+l)*sumint; 

% Calculate right hand side sum over latitude lines if m = 0 
if m == 0 
sumint=0; 
for j=l:t 
integ=0; 

phi=pi*j/(t+l)-eff_thickness/r(b)/2; 

dphi=eff_thickness/r(b)/divs; 
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for i=l:divs 

L=legendre(k,cos(phi) , 'norm') ; 
integ=integ+L( 1 )*gsinc(phi- 
pi*j/(t+l),eff_thickness,r(b),numstddevs)*sin(phi)*dphi; 
%legendre(0,cos(phi),'norm')*g(phi, thickness, r(b))*sin(phi)*dphi; 
phi=phi+dphi; 
end 

sumint=sumint+integ; 

end 

RHS(b,k+l,m+l) = RHS(b,k+l,m+l) + 2*pi*sumint; 
end 

end % close m loop 
%end % close if k is even check 
end % close k loop 
end % close b loop 

for k=0:kmax 
for m=0:k 

% Calculate coefficients: 

Abar(n+ 1 ,k+ 1 ,m+ 1 )=0; 

Bbar( 1 ,k+ 1 ,m+ 1 )=0; 
for i=l:n 
if i>l 

Abar(n+ 1 -i,k+ 1 ,m+ 1 )=Abar(n+ 1 -i+ 1 ,k+ 1 ,m+ 1 )*(r(n+ 1 -i)/r(n+ 1 -i+ 1 )) A (k) + 
q(n+ 1 -i)/thickness/(2*(2*k+ 1 ) *pi*epso)*RHS(n+ 1 -i,k+ 1 ,m+ 1 ) *r(n+ 1 -i); %r(n+ 1 -i) A k 
else 

Abar(n+ 1 -i,k+ 1 ,m+ 1 )= q(n+ 1 -i)/thickness/(2*(2*k+ 1 ) *pi*epso)*RHS(n+ 1 - 
i,k+ 1 ,m+ l)*r(n+ 1 -i); 
end 
if i<n 

Bbar(i+ 1 ,k+ 1 ,m+ 1 )=(Bbar(i,k+ 1 ,m+ 1 ) + 
q(i)/thickness/(2*(2*k+ 1 ) *pi*epso) *RHS(i,k+ 1 ,m+ 1 ) *r(i)) *(r(i)/r(i+l)) A (k+l); 
else 

Bbar(i+ 1 ,k+ 1 ,m+ 1 )=(Bbar(i,k+ 1 ,m+ 1 ) + 
q(i)/thickness/(2*(2*k+ 1 ) *pi*epso)*RHS(i,k+ 1 ,m+ l)*r(i)) ; 
end 
end 
end 
end 

%close(wait); 
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spherical_potential_norm_posm_bar.ni 

% spherical_potential.m 

% This is a Matlab function to calculate the potential at a given position 
% described in spherical coordinates (r, theta, phi) where r is the radius, 

% theta is the longitude angle (0-2*pi), and phi is the co-latitude (O-pi). 

% Written by: Carl Dietrich 8/6/04 

% Inputs: A and B are the coefficient matrices as calculated by the 
% function 'phi_coefficients.m'. The other inputs are described above. 

% Output: This function outputs the scalar electrostatic potential at the 
% specified postition. 

function pot = spherical_potential_norm_posm_bar(rvect,Abar,Bbar,r, theta, phi) 
[nmag,kmag,mmag] = size(Abar); 
n = nmag-1; 
i=l; 

if r>rvect(max(size(rvect))) 
i = max(size(rvect)); 
else 

while r>rvect(i) 
i=i+l; 
end 
end 

pot = 0; 

for k = 0:(kmag-l) 

if k==0 I rem(k,2)==0 %only even ks contribute (assumes n-s symmetry) 

LP = legendre(k,cos(phi),'norm'); 
for m = 0:k 

pot = pot+(Abar(i,k+l,m+l)*(r/rvect(i)) A k+Bbar(i,k+l,m+l)*(r/rvect(i)) A (-k- 
1 )) *LP(m+ l)*cos(m*theta) ; %/((- 1 ) A m*sqrt((k+ 1 )/2*factorial(k-m)/factorial(k+m))) 
end 
end 
end 
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OOPIC Code 

GridDiagnostics5 . inp 

{ 

GridDiagnostics4 . inp — turns off injection at low density and removes 
absorbing injector 

GridDiagnosticsl . inp — based on can.of wormsl 3b, removed: other 1 Injectors 
canofwormsl3b — off axis argon injection at 3.2uA, nonabsorbing 
injector 

canofworms5 mess with grid potentials (4 was shrunk size) 

canofworms3 move emitters in and to lower potential 

canofworms2 lower voltages, for experiment, -Jul y 10, seems to work 

astalkll.inp add a pihi diagnostic 

astalk9.inp add ii 3rd and 4th stalk 

astalk8.inp add in 2nd stalk for other grids 

astalk7 addin a user-defined diagnostic 

astalk6 add in fourth grid 

astalk5 add in a third grid 

astalk4 add in a second grid 

astalk3 go back to ions 

astalk2 add a stalk and fix grid shift problem, also implement 
electron injection 

alt9c add additional anode grids near diagonals 
alt 9b turn, off diagonal injection, add secondary conditions to: 
regions 

alt9 use an 80 by 80 grid with electrodes on angles 
alt6 change this file for each case 

alt5 alternative grid config, begin initial matrix, base case, Aug. 
26, 2003 

ver20 add in dielectric regions inside grid boxes 
verl9 alternate grid structure 

verl5 addi lift other six grids, SHOWS AWESOME FOCUSING ABILITY 
verl4 add in second additional grid 

verl3 expand to a 40x40 grid, traditional IEC version 

ver6 add solid grids 

ver5 add pre-loaded ions 

ver4 add variables 

ver3 adopt quarter config, 20x20 

ver2 adopt el ectror.induced ionization, verify reflections 
verl adapt plasma ignition for cylindrical IEC 

} 

Variables 

{ 

anode = 0 / /anode wall potential in volts 
gridA = -7500 //trapping grid ptoential in volts 
gridB = -10000 //cathode grid potential in volts 
gridC = -5000 //pinching grid near cathode 
gridD = -2000 //pinching grid near anode 
//gridE = 4000 //pinching grid near anode 

curr = 20e-6 //gun currents in Amps 

// idrift = 2e7 //electron injection velocity, to penetrate XOOkeV 
field about 1 cm 

idrift = 3e2 // used for ion injection velocity about 100 eV for D 

seeprob =.8 //I //0 
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wght >= 1 e5 //5e6 //particle weighting 

gdm = 1 //emitter cell size 

dem = 25 //diagonal emitter location. 

oem = 35 //orthogonal emitter location 

dFAC = 0 //turn diagonal injection on or off 

oFAC = 1 //turn orthogonal injection on or off 

pul sedurat i on = le-5 / /2e-3 // ion gun pulse durat 

A1 = 10 //initial x-position of gridA 

A2 = 4 //initial y-positi*>h of gridA 

B1 = 15 

B2 = 6 

Cl = 20 

C2 = 8 

Si = 27 

D2 = 11 

El = 30 

E2 = 12 

FI = 37 //anode grid position 
F2 = 15 

G1 = 34 //additional anode grid segments 
G2 **.' M 

HI = 30 //additional anode grid segments 
H2 = 27 


} 

Region 

{ 

Grid 

{ 

J = 80 
xls = 0,0 
xlf = 0.4 

«a = i . o 

K = 80 
x2s = 0.0 
x2f = 0.4 
n2 = 1.0 
Geometry = 1 

} 

Control 

{ 

dt = 5.0E-9 

F,1 ectrostati cFl ag = 1 

BoltzmannFlag = 0 

} 

MCC 

{ 

gas = Ar //H 
pressure = 3e-8 //3e-6 
eSpecies = electrons 
i Species s=- ions 
collisionFlag = 1 

} 

Species 

{ 
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name - electrons 
m = 9 . HE— 3.1 
q = -1 . 6e-l 9 
particleLimit = 2e6 
collisionModel = 1 
supercycle = 50 

} 

Species 

{ 

name = ions 
m == 6.68e-26 // Argon' 
q = 1.6e-19 // single ionization 
particleLimit = 2e6 
collisionModel = 2 

} 

/// /////////////////////// // ///////////ION 
GUN:///////////////////////////////////////// 

SaitPort //left emitter , vertical 

{ 

speciesName = ions 
'*= 40-oem 

kl = 40 - 2* gdm //37 //40 - gdm 
52 » 40-oem 

k2 = 40 //38 //40 + gdm 
normal = 1 

I = curr*oFAC 

al=l 

a0=0 

trise=0 

tpulse = pulseduration /./ cut off the injection after this much 
tfall=0 

vldrift = idrift 
v2drift = 0 
np2c = wght 
Secondary 
{ 

secondary = secprob 
secSpecies = ions 
iSpecies = ions 

} 

} 

// Now make sidewalls of ion gun: 

Dielectric 

{ 

QuseFlag = 1 
jl = 0 
kl = 40 
j2 = 40 - oem 
kt = 40 
normal = 1 

} 

Dielectric 

{ 

QuseFlag = 1 
jl = 0 

kl = 40 - 2* gdm 
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j2 ; =40- oeir. 
k2 = 40 - 2*gdm 
normal = -1 

} 


/////////////////////////////////////////////// 

Stalk:////////////////////////////////// 

Dielectric //left side of middle STALK 

{ 

QuseFlag = 0 
j$ = 38 
kl = 0 
j2 = 39 
k| = 0 
normal = 1 

} 

Dielectric //right side of middleSTALK 

{ 

QuseFlag = 0 
jl = 41 
kl = 0 
j2 = 42 
k2 = 0 
normal = 1 


Dielectric 

{ 

QuseFlag = 
jl = 36 
kl = 0 
j2 = 37 
k2 = 0 


//l eft side of 2nd STALK 
0 


} 

Dielectric //right side of 2nd STALK 

{ 

QuseFlag = 0 
jl = 43 
kl = 0 
j2 = 44 
k2 = 0 
normal = 1 

} 


Dielectric //left side of 3rd STALK 

{ 

QuseFlag = 0 
jl = 34 
kl = 0 
j2 = 35 
k2 = 0 
normal = 1 

} 

Dielectric //right side of 3rd STALK 

{ 

QuseFlag = 0 
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kl = 80 
j2 = 80 
k2 = 80 
normal = -1 
} 

ft Segment //bottom wall 

// { 

// jl = 0 
// kl; fe 0 
// j2 = 80 
// k2 = 0 
/ / normal = 1 
// } 

Segment //bottom right wa 1 1 //must adjust for stalks 



kl = 0 
j2 = 32 
k2 = 0 



Segment //bottom left wall //must adjust for stalks 



kl = 0 
j2 = 80 
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k2 = 0 
normal = 
} 

Segment 

{ 

jl = 0 
kl = 0 
j2 = 0 
k2 = 80 


//left wall 


} 

Segment //right wall-top 

{ 

jl = 80 

kl = 0 

j2 = 80 

k2 = 80 

normal = -1 

} 

// segments for anode grid 1 
Segment //top wall 
{ 

jl = 40 + FI 

kl = 40 + F2 + 1 

j2 = 40 + FI + 1 

k2 = 40 + F2 + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl =: 40 + FI 

kl = 40 + F2 

j2 = 40 + FI + 1 

k2 = 40 + F2 

normal = -1 

} 

Segment //right wall 

{ 

jl = 40 + FI + 1 
kl = 40 + F2 
j2n*#v4 0 + FI + 1 
k2 = 40 + F2 + 1 
normal = 1 
} 

Segment // left wall 

{ 

jl = 40 + FI 

kl = 40 + F2 

j2 = 40 + FI 

k2 = 40 + F2 + 1 

normal 4 -1 


} 

// segments for anode grid 2 
Segment //top wall 
{ 
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jl = 40 + F2 

kl = 40 + FI + 1 

jg = 40 + F2 + 1 

k2 = 40 + FI + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl = 40 + F2 

kl = 40 + FI 

j2 = 40 + F2 + 1 

k2 = 40 + FI 

normal =. -1 

Segment //right wall 

{ 

jl = 40 + F2 + 1 

kl = 40 + FI 

j2 = 40 + F2 + 1 

k2 = 40 + FI + 1 

normal 1 

Segment // left wall 

{ 

jl = 40 + F2 

kl = 40 + FI 

j2 = 40 + F2 

k2 = 40 + FI + 1 

normal = -1 

} 

// segments for anode grid 3 
Segment //top wall 
{ 

jl = 40 - F2 

kl = 40 + FI + 1 

j2 = 40 - F2 - 1 

k2 = 40 + FI + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl = 40 - F2 

kl = 40 + FI 

j2 = 40 - F2 - 1 

k2 = 40 + FI 

normal = -1 

} 

Segment //right wall 

{ 

jl = 40 - F2 

kl = 40 + FI 

jf = 40 - F2 

k2 =; 40 + FI + 1 

normal = 1 

} 

Segment // left wall 

{ 
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jl = 40 - F2 - 1 
kl = 40 + FI 
jg = 40 - F2 - 1 
k2 = 40 + FI + 1 
normal = -1 
} 

// segments for anode grid 4 
Segment //top wall 
{ 

jl = 40 - FI 
kl = 40 + F2 + 1 
j2 •=• 4 0 - FI - 1 
k2 = 40 + F2 + 1 
normal = 1 
} 

Segment //bottom wall 

{ 

jl = 40 - FI 

kl = 40 + F2 

j2 = 40 - FI - 1 

k2 = 40 + F2 

normal = — 1 

} 

Segment //right wall 

{ 

jl = 40 - FI 

kl = 40 + F2 

j2 = 40 - FI 

k2 = 40 + F2 + 1 

normal = 1 

} 

Segment // left wall 

{ 

jl = 40 - FI - 1 

kl = 40 + F2 

j2 = 40 - FI - 1 

k2 = 40 + F2 + 1 

normal = -1 

} 

// segments for anode grid 8 
Segment //top wall 
{ 

jl = 40 + FI 

kl = 40 - F2 

j2 = 40 + FI + 1 

k2 = 40 - F2 

normal <d 1 

} 

Segment //bottom wall 

{ 

jl = 40 + FI 
kl = 40 - F2 - 1 
j2 =40 • FI + 1 
k2 = 40 - F2 - 1 
normal = -1 
} 

Segment / / right wall 
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{ 

jl = 40 + FI + 1 
kl = 40 - F2 
j2 = 40 + FI + 1 
k2 = 40 - F2 - 1 
normal = 1 
} 

Segment // left wall 

{ 

jl = 40 + FI 

kl = 40 - F2 

j2 = 40 + FI 

k2 if 40 - F2 - 1 

normal = -1 

} 

// segments for anode grid 7 
Segment //top wall 
{ 

jl = 40 + F2 

kl = 40 - FI 

j2 = 40 + F2 + 1 

k2;4-- 4 0 - FI 
normal = 1 
} 

Segment //bottom wall 

{ 

jl = 40 + F2 

kl = 40 - FI - 1 

j2 = 40 + F2 + 1 

k2 = 40 - FI - 1 

normal = -1 

} 

Segment //right wall 

{ 

jl = 40 + F2 + 1 

kl = 40 - FI 

j2 s 40 + F2 + 1 

k2 = 40 - FI - 1 

normal = 1 

} 

Segment // left wall 

{ 

jl = 40 + F2 
kl = 40 - FI 
j2 = 40 + F2 
k2 40 - FI - 1 
normal = -1 
} 

// segments for anode grid 6 
Segment //top wall 
{ 

jl =' 40 - F2 
kl = 40 - FI 
j2 = 40 - F2 - 1 
k2 = 40 - FI 
normal = 1 
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} 

Segment //bottom wall 

{ 

jl = 40 - F2 
kl 40 - FI - 1 
j2 = 40 - F2 - 1 
k2 40 - FI - 1 
normal = -1 
} 

Segment //right wall 

{ 

jl = 40 - F2 

kl =-■ 40 - FI 

j2 = 40 - F2 

k2 = 40 - FI - 1 

normal = X 

} 

Segment // left wall 

{ 

jl = 40 - F2 - 1 
kl = 40 - FI 
j2 = 40 - F2 - 1 
k2;#40 - FI - 1 
normal = -1 
} 

// segments for anode grid 5 
Segment //top wall 
{ 

jl = 40 - FI 
kl = 40 - F2 
j2 = 40 - FI - 1 
k2 =? 40 - F2 
normal = 1 
} 

Segment //bottom wall 

{ 

jl = 40 - FI 

kl = 40 - F2 - 1 

j2 = 40 - FI - 1 

k2 = 40 - F2 - 1 

normal = -1 

} 

Segment //right wall 

{ 

jl = 40 - FI 

kl = 40 - F2 

j2 40 - FI 
k2 m 40 - F2 - 1 
normal = 1 
} 

Segment // left wall 

{ 

jl =' 40 - FI -1 

kl = 40 - F2 

j2 = 40 - FI - 1 

k2 = 40 - F2 - 1 

normal = -1 
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} 

//top right diagonal anode grid segments, 1-4 
// segments for anode grid G, top right lower 
Segment //top wall 
{ 

jl = 40 + G1 

kl /=: 40 + G2 + 1 

j2 = 40 + G1 + 1 

k2 = 40 + G2 + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl = 40 + G1 

kl = 40 + G2 

j2 = 40 + G1 + 1 

k2 *=■ 40 + G2 

normal = -1 

} 

Segment / / right wall 

{ 

jl = 40 + G1 + 1 

kl = 40 + G2 

j2 = 40 + G1 + 1 

k2 = 40 + G2 + 1 

normal = 1 

} 

Segment // left wall 

{ 

jl = 40 + G1 

kl = 40 + G2 

j2 = 40 + G1 

k2 = 40 + G2 + 1 

normal = -1 

} 

// segments for anode grid G, top right upper 
Segment //top wall 
{ 

jl = 40 + G2 

kl = 40 + G1 + 1 

j2 = 40 + G2 + 1 

k2 = 40 + G1 + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl = 40 + G2 

kl = 40 + G1 

j2 = 40 + G2 + 1 

k2 = 40 + G1 

normal = -1 

} 

Segment //right wall 

{ 

jl = 40 + G2 + 1 
kl = 40 + G1 
j2 = 40 + G2 + 1 
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k2 = 40 + G1 + 1 
normal = 1 
} 

Segment // left wall 

{ 

jl = 40 + G2 

kl = 40 + G1 

j2 = 40 + G2 

k2 = 40 + G1 + 1 

normal = -1 

} 

// segments for anode grid H, top right lower 
Segment //top wall 
{ 

jl = 40 + HI 
kl = 40 + ip + 1 
j2 = 40 + HI + 1 
k2 =40 + m + 1 
normal = 1 
} 

Segment //bottom wall 

{ 

jl = 40 + HI 
kl = 40 + H2 
j2 = 40 + HI + 1 
k2 = 40 + H2' 
normal = -1 
} 

Segment //right wall 

{ 

jl = 40 + HI + 1 

kl = 40 + H2 

jt = 40 + HI + 1 

k2 = 40 + H2 + 1 

normal = 1 

} 

Segment // left wall 

{ 

jl = 40 + HI 

kl = 40 + H2 

j'2 = 40 + HI 

k2 = 40 + H2 + 1 

normal = ?-l 

} 

// segments for anode grid H, top right upper 
Segment //top wall 
{ 

jl = 40 + I!2 

kl = 40 + HI + 1 

j2 = 40 + H2 + 1 

k2 = 40 + HI. + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl = 40 + Ip 
kl = 40 + HI 
j2 40 • H2 + 1 
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k2 = 40 + HI 
normal = -1 
} 

Segment //right wall 

{ 

jl = 40 + 12 + 1 
kl = 40 + HI 
j2 = 40 + H2 + 1 
k2 = 40 + HI + 1 
normal = 1 
} 

Segment // left wall 

{ 

jl = 40 + H2 

kl = 40 + HI 

j:2- = 40 + m 

k2 = 40 + HI + 1 

normal = -1 

} 

//top LEFT diagonal anode grid segments, 1-4 
// segments for anode grid G, top LEFT lower 
Segment //top wall 
{ 

jl = 40 - G1 

kl = 40 + G2 + 1 

j2 = 40 - G1 - 1 

k2 = 40 + G2 + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl U 40 - G1 

kl = 40 + G2 

j2 = 40 - G1 - 1 

k2 = 40 + G2 

normal = -1 

} 

Segment //right wall 

{ 

jl = 40 - G1 - 1 
kl = 40 + G2 
j2 '<=;• 4 0 - G1 - 1 
k2 = 40 + G2 + 1 
normal = -1 
} 

Segment // left wall 

{ 

jl = 40 - G1 

kl = 40 + G2 

j2 = 40 - G1 

k2 = 40 + G2 + 1 

normal = 1 

} 

// segments for anode grid G, top LEFT upper 
Segment //top wall 
{ 

jl = 40 - G2 
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kl = 40 + G1 + 1 

j2 = 40 - G2 - 1 

k2 = 40 + G1 + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl m 40 - G2 

kl = 40 + G1 

j2 = 40 - G2 - 1 

k2 = 40 + G1 

normal = -1 

} 

Segment //right wall 

{ 

jl = 40 - G2 - 1 

kl = 40 + G1 

j2 = 40 - G2 - 1 

k2 = 40 + G1 + 1 

normal = -1 

} 

Segment // left wall 

{ 

jl = 40 - G2 

kl = 40 + G1 

j2 = 40 - G2 

k2 = 40 + G1 + 1 

normal = 1 

} 

// segments for anode grid H, top LEFT lower 
Segment //top wall 
{ 

jl =■ 40 - HI 
kl = 40 + H2 + 1 

j2 = 40 - HI - 1 

k2 40 • H:2 + 1 
normal = 1 
} 

Segment //bottom wall 

{ 

jl = 40 - HI 

kl = 40 + H2 

j2 '<=;• 4 0 - HI - 1 

k2 = 40 + H2 
normal = -1 
} 

Segment //right wall 

{ 

jl = 40 - HI - 1 

kl = 40 + H2 

j2 = 4 0 - HI. - 1 

k2 = 40 + H2 + 1 

normal = -1 

} 

Segment // left wall 

{ 

jl = 40 - HI 

kl = 40 + H2 
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j2 = 40 - HI 
k2 = 40 + H2 + 1 
normal = 1 
} 

// segments for anode grid H, top LEFT upper 
Segment //top wall 
{ 

jl = 40 - H2 

kl = 40 + HI + 1 

j2 = 40 - H2 - 1 

k2 = 40 + HI. + 1 

normal = 1 

} 

Segment //bottom wall 

{ 

jl = 40 - ,1:2; 

kl = 40 + HI 

j2 = 40 - 12. - 1 

k2 = 40 + HI 

normal = -1 

} 

Segment //right wall 

{ 

jl = 40 - H2 - 1 
kl = 40 + HI 
j2 = 40 - ,12 - 1 
k2 = 40 + HI + 1 
normal = -1 
} 

Segment // left wall 

{ 

jl = 40 - H2 
kl if 40 + HI 
j2 = 40 - H2 
k2 = 40 + HI + 1 
normal = 1 
} 


// segments for anode grid G, BOTTOM right lower 
Segment //top wall 
{ 

jl = 40 + G1 

kl = 40 - G2 - 1 

j2 = 40 + G1 + 1 

k2 = 40 - G2 - 1 

normal = -1 

} 

Segment //bottom wall 

{ 

jl = 40 + G1 

kl = 40 - G2 

ji = 40 + G1 + 1 

k2 =;• 4 0 - G2 

normal = 1 

} 

Segment //right wall 

{ 
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jl = 40 + G1 + 1 
kl = 40 - G2 
j'£ = 40 + G1 + 1 
k2 = 40 - G2 - 1 
normal = 1 
} 

Segment // left wall 

{ 

jl = 40 + G1 

kl = 40 - G2 

j2 = 40 + G1 

k2 = 40 - G2 - 1 

normal = *-1 

} 

// segments for anode grid G, BOTTOM right upper 
Segment //top wall 
{ 

jl = 40 + G2 

kl = 40 - G1 - 1 

j2 = 40 + G2 + 1 

k2 = 40 - G1 - 1 

normal = -1 

} 

Segment //bottom wall 

{ 

jl = 40 + G2 

kl = 40 - G1 

j2 = 40 + G2 + 1 

k2 = 40 - G1 

normal = 1 

} 

Segment //right wall 

{ 

jl = 40 + G2 + 1 

kl = 40 - G1 

j2 = 40 + G2 + 1 

k2 = 40 - G1 - 1 

normal == 1 

} 

Segment // left wall 

{ 

jl = 40 + G2 

kl = 40 - G1 

j2 = 40 + G2 

k2 = 40 - G1 - 1 

normal = -1 

} 

// segments for anode grid H, BOTTOM right lower 
Segment //top wall 

{ 

jl = 40 + HI 
kl = 40 - H2 - 1 
jl' = 4 0 + HI + 1 

k2 =;• 4 0 - H2 - 1 

normal = -1 

} 

Segment //bottom wall 

{ 
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jl = 40 + HI 
kl = 40 - H2 
j'£ = 40 + HI + 1 
k2 = 40 - H2 
normal = 1 
} 

Segment //right wall 

{ 

jl = 40 + HI + 1 

kl = 40 - H2 

■j 2 = 40 + HI. + 1 

k2 = 40 - H2 - 1 

normal = 1 

} 

Segment // left wall 

{ 

jl = 40 + HI 

kl = 40 - H2 

j2 = 40 + HI 

k2 = 40 - H2 - 1 

normal = -1 

} 

// segments for anode grid H, BOTTOM right upper 
Segment //top wall 
{ 

jl = 40 + 0 
kl = 40 - HI - 1 
j2 4 0 + fH + 1 

k2 = 40 - HI - 1 
normal = -1 
} 

Segment //bottom wall 

{ 

jl = 40 + H2 
kl = 40 - HI 
j2 4 0 • H2: + 1 

k2 = 40 - HI 
normal ~ 1 
} 

Segment //right wall 

{ 

jl = 40 + H2 + 1 

kl = 40 - HI 

j2 = 40 + H2 + 1 

k2 = 40 - HI - 1 

normal = 1 

} 

Segment // left wall 

{ 

jl = 40 + H2 

kl = 40 - HI 

j2 = 40 + H2 

k2 '.te. 40 - HI - 1 

normal = -1 

} 
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//BOTTOM LEFT diagonal anode grid segments, 1-4 
// segments for anode grid G, BOTTOM LEFT lower 
Segment //top wall 
{ 

jl = 40 - G1 
kl = 40 - G2 -. r 1 
j2 = 40 - G1 - 1 
k2.' = 40 - G2 - 1 
normal = -1 
} 

Segment //bottom wall. 

{ 

jl = 40 - G1 
kl = 40 - G2 
j2 = 40 - G1 - 1 
k2 = 40 - G2 
normal = 1 
} 

Segment //right wall 

{ 

jl * 40 - G1 - 1 

kl = 40 - G2 

j 2 = 40 - G1 - 1 

k2 = 40 - G2 - 1 

normal = -1 

} 

Segment // left wall 

{ 

jl = 40 - G1 

kl = 40 - G2 

j2- = 40 - G1 

k2 = 40 - G2 - 1 

normal 1 

} 

// segments for anode grid G, BOTTOM LEFT upper 
Segment //top wall 
{ 

jl = 40 - G2 

kl = 40 - G1 - 1 

j2 = 40 - G2 - 1 

kf: i 40 - G1 - 1 

normal = -1 

} 

Segment //bottom wall 

{ 

jl = 40 - G2 

kl = 40 - G1 

j2 == 40 - G2: - 1 

kt « 40 - G1 

normal = 1 

} 

Segment //right wall 

{ 

jl =■ 40 - G2 ~ 1 

kl = 40 - G1 

j2 = 40 - G2 - 1 

k2 = 40 - G1 - 1 

normal = -1 
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} 

Segment // left wall 

{ 

jl = 40 - G2 

kl = 40 - G1 

j2 = 40 - G2 

k2 = 40 - G1 - 1 

normal = 1 

} 

// segments for anode grid H, BOTTOM LEFT lower 
Segment //top wall 
{ 

jl = 40 - HI 

kl = 40 - H2 - 1 

j2 = 40 - HI - 1 

k2 40 - .HZ - 1 
normal = -1 
} 

Segment //bottom wall 

{ 

jl = 40 - HI 

kl = 40 - H2 

j2 -4,40 - HI - 1 
k2 = 40 - H2 
normal = 1 
} 

Segment //right wall 

{ 

jl = 40 - HI - 1 
kl = 40 - H2 
j2- = 40 - HI - 1 
k2 = 40 - H2 - 1 
normal = -1 
} 

Segment // left wall 

{ 

jl = 40 - HI 

kl = 40 - HZ 

j2 = 40 - HI 

k2 = 40 - H2 - 1 

normal % % 

} 

/ / segments for anode grid H, BOTTOM LEFT upper 
Segment //top wall 
{ 

jl = 40 - HZ 
kl = 40 - HI - 1 
j.2 m 4 0 - HZ: - 1 
k2 m 40 - HI - 1 
normal = -1 
} 

Segment //bottom wall 

{ 

jl =' 40 - H2 
kl = 40 - HI 
jZ = 40 - HZ - 1 
k2 = 40 - HI 
normal = 1 
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} 

Segment //right wall 

{ 

jl = 40 - H2 - 1 
kl = 40 - HI 
j2 4 0 - - 1 

k2 = 40 - HI - 1 
normal = -1 
} 

Segment // left wall 

{ 

jl = 40 - H2 

kl =-■ 40 - HI 

j2 = 40 - H2 

k2 = 40 - HI - 1 

normal = 1 

} 


////////////////////////////////// End of ANODE* Start Grid 
A //////////////////////// 


Equipotential //gridA 

{ 

C = gridA 
Secondary 
{ 

secondary = 1 
secSpecies = electrons 
iSpecies = ions 

} 

//grid 1 

Segment //top 1 

{ 

jl = 40 + A1 

kl = 40 + A2 + 1 

j2 = 40 + A1 + 1 

k2 = 40 + A2 + 1 

normal = 1 

} 

Segment //bottom 1 

{ 

jl = 40 + A1 

kl = 40 + A2 

j2 = 40 + A1 + 1 

k2 = 40 + A2 

normal = -1 

} 

Segment //right 1 

{ 

jl = 40 + A1 + 1 
kl = 40 + A2 
j2 = 40 + A1 + 1 
k2- = 40 + A2 + 1 
normal = 1 
} 
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Segment // left 1 

{ 

jl = 40 + A1 

kl = 40 + A2 

j2 = 40 + A1 

k2 = 40 + A2 + 1 



//grid 2 

Segment //top 1 

{ 

jl = 40 + A2 

kl ** 40 + A1 + 1 

j2 = 40 + A2 + 1 

k2 = 40 + A1 + 1 

normal *. .1 

Segment //bottom 1 

{ 

jl = 40 + A2 

kl = 40 + A1 

j2 = 40 + A2 + 1 

k2 = 40 + A1 



% 40 • A2 + 

k2 = 40 + A1 + 



Segment // left 1 

{ 

I' = 40 + A2 



//grid 4 

Segment / /top 1 

{ 


jl = 40 - A1 

kl = 40 + A2 + 1 

j.2 40 - a: 

k2 = 40 + A2 + 1 

norma 1 = 1 

} 

Segment //bottom 1 

{ 

W = 40 - A1 

kl = 40 + A2 

j-2 = 40 - A1 - 1 

k2 = 40 + A2 

normal = -1 

} 
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Segment / / right 1 

{ 

= 4 0 - A1 

kl = 40 + A2 

j?. 4 0 a: 

k2 = 40 + A2 + 1 

normal = 1 

} 

Segment // left 1 

{ 

rM = 40 - A1 - 1 

kl = 40 + A2 

,j2 = 40 - A1 - 1 

k2 = 40 + A2 + 1 

normal = -1 

} 

//grid 3 

Segment //top 1 

{ 

jl = 40 - A2 

kl = 40 + A1 + 1 

j2 = 40 - A2 - 1 

k2 = 40 + A1 + 1 

normal = 1 

} 

Segment //bottom 1 

{ 

= 4 0 - A2 
kl = 40 + A1 
j?. = 40 - A2 - 1 
k2 = 40 + A1 
normal = -1 
} 

Segment //right 1 

{ 

4,4; = 4 0 - A2 
kl = 40 + A1 
;J2' = 4 0 - A2 
k2 = 40 + A1 + 1 
normal = 1 
) 

Segment // left 1 

{ 

U = 40 - A2 -1 

kl = 40 + A1 

= 40 - A2 -1 
k2 = 40 + A1 + 1 
norma 1 = -1 
} 

//grid 8 

Segment / /top 1 

{ 

ia = 40 + A1 

kl = 40 - A2 

j2 = 40 + A1 + 1 

k2 = 40 - A2 

normal = 1 

} 
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Segment //bottom 1 

{ 

jf?L = 40 + A1 

kl = 40 - A2 - 1 

j2 = 40 + A1 + 1 

k2 = 40 - A2 - 1 

normal = -1 

} 

Segment //right 1 

{ 

;|a = 40 + ai + i 

kl = 40 - A2 

jf2 = 40 + Al + 1 

k2 = 40 - A2 - 1 

normal = 1 

} 

Segment // left 1 

{ 

jl = 40 + Al 

kl = 40 - A2 

!|2 = 40 + Al 

k2 = 40 - A2 - 1 

norma 1 - -1 

} 

//grid 7 

Segment //top 1 

{ 

;Jt = 40 
kl = 40 - Al 
j2 = 40 + A2 + 1 
k2 = 40 - Al 
normal = 1 
} 

Segment //bottom 1 

{ 

l|a = 41 

kl = 40 - Al - 1 

;J2 = 40 + A2 + 1 

k2 = 40 - Al - 1 

normal = -1 

4 

Segment //right 1 

{ 

jl = 40 + A2 + 1 

kl = 40 - Al 

= 40 + A2 + 1 

k2 = 40 - Al - 1 

normal = 1 

} 

// Segment // left 1 

fj { 

// jl = 40 + A2 

// kl = 40 - Al 

// = 40 + A2 

// k2 = 40 - Al + 1 
// normal = -1 

// } 

//grid 5 
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Segment / /top 1 

{ 

yif = 40 - A1 

kl = 40 - A2 

j2 = 40 - A1 - 1 
k2 = 40 - A2 
normal = 1 
} 

Segment //bottom 1 

{ 

;|a = 40 - ai 

kl = 40 - A2 -1 

. m = 40 - A1 - 1 

k2 = 40 - A2 -1 

normal = -1 

} 

Segment //right 1 

{ 

jl = 40 - Al 

kl = 40 - A2 

!|2 = 40 - Al 

k2 = 40 - A2 - 1 

norma 1 = 1 

} 

Segment // left 1 

{ 

jl = 40 - Al -1 

kl = 40 - A2 

j2 = 40 - Al -1 

k2 = 40 - A2 - 1 

normal = -1 

} 

//grid 6 

Segment //top 1 

{ 

= 40 

kl = 40 - Al 
;$2 = 40 - A2 - 1 
k2 = 40 - Al 
normal = 1 
} 

Segment //bottom 1 

{ 

jl = 39 

kl = 40 - Al - 1 

j|2 = 40 - A2 - 1 

k2 = 40 - Al - 1 

normal = -1 

} 

// Segment //right 1 

fj { 

// jl = 40 - A2 + 1 

// kl = 40 - Al 

// ^ = 40 - A2 + 1 

// k2 = 40 - Al + 1 

// normal = 1 

// } 

Segment // left 1 
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{ 

jl = 40 - A2 - 1 

kl = 40 - A1 

j2 = 40 - A2 - 1 

k2 = 40 - A1 - 1 

normal = -1 

i 

//add stalk region 

Segment // left stalk 

{ 

3§§ = 39 

kl = 0 

p = 39 

k2 = 40 - A1 - 1 

normal = -1 

} 

Segment // right stalk 

{ 

jl = 41 
kl = 0 
j2 = 41 

k2 = 40 - A1 - 1 
normal = 1 
} 

Segment // bottom stalk 

{ 

Ip = 39 

kl = 0 
J|| = 41 
k2 = 0 
normal = 1 
} 

} 

//////////////////////////////////End of grid A: Start Grid 
B ://///////////////////// 

// grid B 

Equipotential //gridB 

{ 

C = gridB 
Secondary 
{ 

secondary = t 
secSpecies = electrons 
i Species = ions 

} 

//grid 1 

Segment / /top 1 

{ 

4l = 40 + B1 
kl = 40 + B2 + 1 
j2 = 40 + B1 + 1 


191 



k2 = 40 + B2 + 1 
normal = 1 
} 

Segment //bottom 1 

{ 

,.|t: = 40 + B1 
kl = 40 + B2 
;|i = 40 + B1 + 1 
k2 = 40 + B2 
normal = -1 
} 

Segment //right 1 

{ 

jl = 40 + B1 + 1 

kl = 40 + B2 

-|-2 = 40 + B1 + 1 

k2 = 40 + B2 + 1 

norma 1 = 1 

} 

Segment // left 1 

{ 

jl = 40 + B1 

kl = 40 + B2 

j-2 40 • 31 

k2 = 40 + B2 + 1 

normal = -1 

} 

//grid 2 

Segment //top 1 

{ 

■;J£ : = 4 0 + B2 
kl = 40 + B1 + 1 
|| = 40 + B2 + 1 
k2 = 40 + B1 + 1 
normal = 1 
} 

Segment //bottom 1 

{ 

jl = 40 + B2 

kl = 40 + B1 

= 40 + B2 + 1 

k2 = 40 + B1 

normal = -1 

} 

Segment //right 1 

{ 

jl = 40 + B2 + 1 

kl = 40 + B1 

= 40 + B2 + 1 

k2 = 40 + B1 + 1 

normal = 1 

| 

Segment // left 1 

{ 

jl = 40 + B2 

kl = 40 + B1 

j;2 = 40 + B2 

k2 = 40 + B1 + 1 
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normal = -1 

} 

//grid 4 

Segment //top 1 

{ 

,.|t; = 40 - B1 

kl = 40 + B2 + 1 

;|i = 4 0 - B1 - 1 

k2 = 40 + B2 + 1 

normal = 1 

} 

Segment //bottom 1 

{ 

jl = 40 - B1 

kl = 40 + B2 

-|-2 = 40 - B1 - 1 

k2 = 40 + B2 

norma 1 = -1 

} 

Segment //right 1 

{ 

jl = 40 - B1 

kl = 40 + B2 

j-2 40 31 

k2 = 40 + B2 + 1 

normal = 1 

} 

Segment // left 1 

{ 

jl = 40 - B1 - 1 

kl = 40 + B2 

j2 = 40 - B1 - 1 

k2 = 40 + B2 + 1 

normal = -1 

} 

//grid 3 

Segment //top 1 

{ 

jl = 40 - B2 
kl = 40 + B1 + 1 
= 40 - B2 - 1 
k2 = 40 + B1 + 1 
normal = 1 
} 

Segment //bottom 1 

{ 

jl = 40 - B2 

kl = 40 + B1 

= 40 - B2 - 1 

k2 = 40 + B1 

normal = -1 

} 

Segment / / right 1 

{ 

jl = 40 - B2 

kl = 40 + B1 

j;2 = 40 - B2 

k2 = 40 + B1 + 1 
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normal = 1 

} 

Segment // left 1 

{ 

jl = 40 - B2 -1 

kl = 40 + B1 

j;2 = 40 - B2 -1 

k2 = 40 + B1 + 1 

normal = -1 

} 

//grid 8 

Segment //top 1 

{ 

jl = 40 + B1 

kl = 40 - B2 

-|-2 = 40 + B1 + 1 

k2 = 40 - B2 

norma 1 = 1 

} 

Segment //bottom 1 

{ 

jl = 40 + B1 

kl = 40 - B2 - 1 

j'2 = 40 + B1 + 1 

k2 = 40 - B2 - 1 

normal = -1 

} 

Segment / / right 1 

{ 

jl = 40 + B1 + 1 

kl = 40 - B2 

j2 = 40 + B1 + 1 

k2 = 40 - B2 - 1 

normal = 1 

} 

Segment // left 1 

{ 

= 40 + B1 

kl = 40 - B2 

j2 40 • 31 

k2 = 40 - B2 - 1 

normal = -1 

1 

//grid 7 

Segment //top 1 

{ 

jl = 40 + B2 

kl = 40 - B1 

|§ = 40 + B2 + 1 

k2 = 40 - B1 

normal = 1 

| 

Segment //bottom 1 

{ 

jl = 40 + B2 

kl = 40 - B1 - 1 

j;2 = 40 + B2 + 1 

k2 = 40 - B1 - 1 
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normal = -1 

} 

Segment //right 1 

{ 

jl = 40 + B2 + 1 

kl = 40 - B1 

j;2 = 40 + B2 + 1 

k2 = 40 - B1 - 1 

normal = 1 

} 

// Segment // left 1 

// { 

// : ||L = 40 + B2 

// kl = 40 - B1 

// J2 = 40 + B2 

// k2 = 40 - B1 - 1 

// normal = -1 

// } 

//grid 5 

Segment //top 1 

{ 

jl = 40 - B1 

kl = 40 - B2 

j2 40-31 
k2 = 40 - B2 
normal = 1 
} 

Segment / /bottom 

{ 

jl = 40 - B1 

kl = 40 - B2 -1 

j2 = 40 - B1 - 1 

k2 = 40 - B2 -1 

normal = -1 

} 

Segment //right 1 

{ 

|| = 40 - B1 

kl = 40 - B2 

j2 40 - 31 
k2 = 40 - B2 - 1 
normal = 1 
1 

Segment // left 1 

{ 

Ja = 40 - B1 -1 

kl = 40 - B2 

j.2 = 40 - B1 -1 

k2 = 40 - B2 - 1 

normal = -1 

} 

//grid 6 

Segment / /top 1 

{ 

jl = 40 - B2 

kl = 40 - B1 

j;2 = 40 - B2 - 1 

k2 = 40 - B1 


//comment out for stalk 
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normal = 1 

} 

Segment //bottom 1 

{ 

jl = 40 - B2 

kl = 40 - B1 - 1 

j;2 = 40 - B2 - 1 

k2 = 40 - B1 - 1 

normal = -1 

} 

// Segment //'right 1 

// { 

// : ||L = 40 - B2 

// kl = 40 - B1 

// J2 = 40 - B2 

// k2 = 40 - B1 - 1 

// normal = 1 

// } 

Segment // left 1 

{ 

;|jL = 40 - B2 - 1 

kl = 40 - B1 

= 40 - B2 - 1 

k2 = 40 - B1 - 1 

normal = -1 

I 

//add B-stalk regions 

Segment // right top 

{ 

4 2 

kl = 40 - B1 

|j£ = 40 + B2 

k2 = 40 - B1 

normal = 1 

} 

Segment // right botto: 

{ 

jl = 43 

kl = 40 - B1 - 1 
;|2 = 40 + B2 
k2 = 40 - B1 - 1 
normal = -1 
} 

Segment // right right 

{ 

jl = 43 
kl = 0 
j2 = 43 

k2 = 40 - B1 - 1 
.normal = 1 
| 

Segment // right left 

{ 

jl = 42 
kl = 0 
j;2 = 42 
k2 = 40 - B1 


//comment out for stalk 
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normal = -1 

} 

Segment / / right bottom 

{ 

jl = 42 
kl = 0 
j:2 = 43 
k2 = 0 
normal = 1 
} 

Segment // left top 

{ 

0 = 40 - B2 
kl = 40 - B1 
j2 = 38 
k2 = 40 - B1 
normal = 1 
} 

Segment // left bottom 

{ 

?§§ = 40 - B2 
kl = 40 - B1 - 1 
}2. = 37 

k2 = 40 - B1 - 1 
normal = -1 
} 

Segment // left right 

{ 

jl = 38 
kl = 0 
j2 = 38 
k2 = 40 - B1 
normal = 1 
} 

Segment // left left 

{ 

jl = 37 
kl = 0 
j2 = 37 

k2 = 40 - B1 - 1 
normal = -1 
} 

Segment // left bottom 

{ 

jl = 37 
kl = 0 
j:2 = 38 
k2 = 0 
normal = 1 
} 


} 

//end of grid B 

///////////////////////////////////End of grid B: Start grid 

C :////////////////////// 
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// grid C 

Equipotential //gridC 

{ 

C = gridC 
Secondary 
{ 

secondary = 1 
secSpecies = electrons 
iSpecies = ions 

} 

//grid 1 

Segment //top 1 

{ 

;Ja = 40 + Cl 

kl = 40 + C2 + 1 

j2 = 40 + Cl + 1 

k2 = 40 + C2 + 1 

normal = 1 

} 

Segment //bottom 1 

{ 

jl = 40 + Cl 

kl = 40 + C2 

i% = 40 + Cl + 1 

k2 = 40 + C2 

normal = -1 

} 

Segment //right 1 

{ 

jl = 40 + Cl + 1 

kl = 40 + C2 

j2 = 40 + Cl + 1 

k2 = 40 + C2 + 1 

normal = 1 

} 

Segment // left 1 

{ 

jl = 40 + Cl 

kl = 40 + C2 

j2 = 40 + Cl 

k2 = 40 + C2 + 1 

normal = -1 

} 

//grid 2 

Segment //top 1 

{ 

*|| = 40 + C2 

kl = 40 + Cl + 1 

:|2 = 40 + C2 + 1 

k2 = 40 + Cl + 1 

normal = 1 

} 

Segment //bottom 1 

{ 

jl = 40 + C2 

kl = 40 + Cl 
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= 40 + C2 + 1 
k2 = 40 + Cl 
normal = -1 
} 

Segment //right 1 


jl = 40 + C2 + 1 



//grid 3 
Segment //top 


jl = 40 - C2 
kl = 40 + Cl 
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k2 = 40 - C2 - 1 
normal = 1 
} 

Segment // left 1 

{ 

JtL = 40 - Cl -1 

kl = 40 - C2 

-j'i> = 40 - Cl -1 

k2 = 40 - C2 - 1 

normal = -1 

) 

//grid 6 

Segment / /top 1 

{ 

jl = 40 - C2 

kl = 40 - Cl 

j:2 = 40 - C2 - 1 

k2 = 40 - Cl 

normal = 1 

} 

Segment //bottom 1 

{ 

= 40 - C2 
kl = 40 - Cl - 1 

j2 = 40 - C2 - 1 

k2 = 40 - Cl - 1 

normal = -1 

| 

// Segment //right 1 

// { 

// = 40 - C2 

// kl = 40 - Cl 

// j2 = 40 - C2 

// k2 = 40 - Cl - 1 

// normal = 1 

// } 

Segment // left 1 

{ 

jl = 40 - C2 - 1 

kl = 40 - Cl 

;|2 = 40 - C2 - 1 

k2 = 40 - Cl - 1 

normal = -1 

} 

//add C-stalk regions 

Segment // right top 

{ 

* = 44 

kl = 40 - Cl 

;j:2 = 40 + C2 

k2 = 40 - Cl 

normal = 1 

} 

Segment // right bottom 

{ 

jl = 45 

kl = 40 - Cl - 1 
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J;2 = 4 0 + C2 
k2 = 40 - Cl - 1 
normal = -1 
} 

Segment // right right 

{ 

jl = 45 
kl = 0 
j2 = 45 

k2 = 40 - Cl - 1 
normal = 1 
} 

Segment // right left 

{ 

jl = 44 
kl = 0 
j2 = 44 
k2 = 40 - Cl 
normal = -1 
} 

Segment / /' right bottom 

{ 

iff =44 
kl = 0 
j2 =45 
k2 = 0 
normal = 1 
} 

Segment // left top 

{ 

; |i4 = 40 - C2 

kl = 40 - Cl 

j2- =36 
k2 = 40 - Cl 
normal = 1 
} 

Segment // left bottom 

{ 

jl 40 - C2 

kl = 40 - Cl - 1 

32 = 35 

k2 = 40 - Cl - 1 

normal = -1 

} 

Segment // left right 

{ 

jl = 36 
kl = 0 
j2 = 36 
k2 = 40 - Cl 
normal = 1 
| 

Segment // left left 

{ 

jl = 35 
kl = 0 
|§ =35 

k2 = 40 - Cl - 1 
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normal = -1 

} 

Segment // left bottom 

{ 

jl = 35 
kl = 0 
j:2 = 36 
k2 = 0 
normal = 1 
} 


} 

//end of grid C 

///////////////////////////////////////End of grid C: Start G§| 

D :////////////////////// 

// grid D 

Equipotential //gridD 

{ 

C = gridD 
Secondary 
{ 

secondary '= 1 
secSpecies = electrons 
dSpecies = ions 

} 

//grid 1 

Segment //top 1 

{ 

jl = 40 + D1 
kl = 40 + D2 + 1 
j2 = 40 + D1 + 1 

k2 = 40 + D2 + 1 

normal = 1 

} 

Segment //bottom 1 

{ 

'M = 40 + D1 

kl = 40 + D2 

; J2 = 40 + D1 + 1 

k2 = 40 + D2 

normal = -1 

} 

Segment //right 1 

{ 

tli = 40 + D1 + 1 

kl = 40 + D2 

;|2 = 40 + D1 + 1 

k2 = 40 + D2 + 1 

normal = 1 

} 

Segment // left 1 

{ 

jl = 40 + D1 

kl = 40 + D2 
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IP = 40 + D1 
k2 = 40 + D2 + 1 
normal = -1 
} 

//grid 2 

Segment / /top 1 

{ 

= 40 + D2 
kl = 40 + D1 + 1 
j& = 40 + D2 + 1 
k2 = 40 + D1 + 1 
normal = 1 
} 

Segment //bottom 1 

{ 

= 40 + D2 

kl = 40 + D1 

jt = 40 + D2 + 1 

k2 = 40 + D1 

normal = -1 

} 

Segment //right 1 

{ 

jl = 40 + D2 + 1 

kl = 40 + D1 

3-2 = 40 + D2 + 1 

k2 = 40 + D1 + 1 

normal = 1 

} 

Segment // left 1 

{ 

jl = 40 + D2 

kl = 40 + D1 

•%% = 40 + D2 

k2 = 40 + D1 + 1 

normal = -1 

} 

//grid 4 

Segment //top 1 

{ 

;g;l' = 4 0 - D1 
kl = 40 + D2 + 1 
= 40 - D1 - 1 
k2 = 40 + D2 + 1 
normal = 1 
} 

Segment //bottom 1 

{ 

'1% = 40 - D1 

kl = 40 + D2 

j2 = 40 - D1 - 1 

k2 = 40 + D2 

normal = -1 

} 

Segment //right 1 

{ 

jl = 40 - D1 

kl = 40 + D2 
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= 4 0 - D1 
k2 = 40 + D2 + 1 
normal = 1 
} 

Segment // left 1 

{ 

jl = 40 - D1 - 1 
kl = 40 + D2 
j?. = 40 - D1 - 1 
k2 = 40 + D2 + 1 
normal = -1 
} 

//grid 3 

Segment //top 1 

{ 

«$&• = 40 - D2 
kl = 40 + D1 + 1 
jt = 40 - D2 - 1 
k2 = 40 + D1 + 1 
normal = 1 
} 

Segment //bottom 1 

{ 

jl = 40 - D2 
kl = 40 + D1 
j2 = 40 - D2 - 1 
k2 = 40 + D1 
normal = -1 
} 

Segment //right 1 

{ 

jl = 40 - D2 

kl = 40 + D1 

j2 = 40 - D2 

k2 = 40 + D1 + 1 

normal = 1 

} 

Segment // left 1 

{ 

jl = 40 - D2 -1 

kl = 40 + D1 

j2 = 40 - D2 -1 

k2 = 40 + D1 + 1 

normal = -1 

} 

//grid 8 

Segment //top 1 

{ 

Ifl; = 40 + D1 

kl = 40 - D2 

j2 = 40 + D1 + 1 

k2 = 40 - D2 

normal = 1 

} 

Segment //bottom 1 

{ 

jl = 40 + D1 
kl = 40 - D2 - 1 
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= 40 + di + l 
k2 = 40 - D2 - 1 
normal = -1 
} 

Segment //right 1 

{ 

jl = 40 + Dl + 1 
kl = 40 - D2 
j?. = 40 + Dl + 1 
k2 = 40 - D2 - 1 
normal = 1 
} 

Segment // left 1 

{ 

jl = 40 + Dl 

kl = 40 - D2 

Jl = 40 + Dl 

k2 = 40 - D2 - 1 

normal = -1 

} 

//grid 7 

Segment //top 1 

{ 

jl = 40 + D2 

kl = 40 - Dl 

j% = 40 + D2 + 1 

k2 = 40 - Dl 

normal = 1 

} 

Segment //bottom 1 

{ 

jl = 40 + D2 

kl = 40 - Dl - 1 

f2 = 40 + D2 + 1 

k2 = 40 - Dl - 1 

normal = -1 

} 

Segment / / right 1 

{ 

jl = 40 + D2 + 1 

kl = 40 - Dl 

j2 = 40 + D2 + 1 

k2 = 40 - Dl - 1 

normal = 1 

} 

Segment // left 1 

{ 

|||, = 40 + D2 

kl = 40 - Dl 

j2 = 40 + D2 

k2 = 40 - Dl - 1 

normal = -1 

} 

//grid 5 

Segment //top 1 

{ 

jl = 40 - Dl 

kl = 40 - D2 
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;|S: = 40 - di - l 
k2 = 40 - D2 
normal = 1 
} 

Segment //bottom 1 

{ 

jl = 40 - Dl 

kl = 40 - D2 -1 

j?. = 40 - Dl - 1 

k2 = 40 - D2 -1 

normal = -1 

} 

Segment / / right 1 

{ 

jl = 40 - Dl 

kl = 40 - D2 

j;2 = 40 - Dl 

k2 = 40 - D2 - 1 

normal = 1 

} 

Segment // left 1 

{ 

|p§. = 40 - Dl -1 

kl = 40 - D2 

J§ = 40 - Dl -1 

k2 = 40 - D2 - 1 

normal = -1 

} 

//grid 6 

Segment //top 1 

{ 

jl = 40 - D2 

kl = 40 - Dl 

j2 = 40 - D2 - 1 

k2 = 40 - Dl 

normal = 1 

} 

Segment //bottom 1 

{ 

jl = 40 - D2 

kl = 40 - Dl - 1 

j;2 = 40 - D2 - 1 

k2 = 40 - Dl - 1 

normal = -1 

} 

Segment / / right 1 

{ 

m = 40 - D2 

kl = 40 - Dl 

j2 = 40 - D2 

k2 = 40 - Dl - 1 

normal = 1 

} 

Segment // left 1 

{ 

Jp; = 40 - D2 - 1 

kl = 40 - Dl 

= 40 - D2 - 1 
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k2 = 40 - D1 - 1 
normal = -1 
} 

//add D-stalk regions 

Segment // right top 

{ 

•Jr! = 46 

kl = 40 - D1 

j2 = 40 + D2 

k2 = 40 - D1 

normal = 1 

} 

Segment // right bottom 

{ 

m = 47 

kl = 40 - D1 - 1 
j2 = 40 + D2 
k2 = 40 - D1 - 1 
normal = -1 
} 

Segment // right right 

{ 

jl = 47 
kl = 0 
jl2 = 47 

k2 = 40 - D1 - 1 
normal = 1 
} 

Segment // right left 

{ 

jl = 46 
kl = 0 
j2 = 46 
k2 = 40 - D1 
normal = -1 
} 

Segment // right bottom 

{ 

jl = 46 
kl = 0 
j2 = 47 
k2 = 0 
normal = 1 
} 

Segment // left top 

{ 

= 40 - D2 
kl = 40 - D1 
j 2 = 34 
k2 = 40 - D1 
normal = 1 
} 

Segment // left bottom 

{ 

4l = 40 - D2 
kl = 40 - D1 - 1 
= 33 
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k2 = 40 - D1 - 1 
normal = -1 
} 

Segment // left right 

{ 

= 34 
kl = 0 
j2 = 34 
k2 = 40 - D1 
normal = 1 
} 

Segment // left left 

{ 

P = 33 
kl = 0 
J2 = 33 

k2 = 40 - D1 - 1 
norma 1 = -1 
} 

Segment // left bottom 

{ 

jl = 33 
kl = 0 
H =34 
k2 = 0 
normal = 1 
} 


//end of grid D 

////////////////////////////////////End of Grid D: Start 'Metric 
Regions////////////// 


//regions for grid A 

DielectricRegion //gridA-1 

{ 

QuseFlag = 0 
jl = 40 + A1 
kl = 40 + A2+1 
j'2 = 40 + A1 
k2 = 40 + A2+1 

} 

DielectricRegion //gridA-2 

{ 

QuseFlag = 0 
jl = 40 + A2 
kl = 40 + Al+1 
j2 = 40 + A2 
k2 = 40 + Al+1 

} 

DielectricRegion //gridA-4 

{ 

QuseFlag = 0 
jl = 40 - A1 
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kl = 40 + A2+1 
j;2 = 40 - A1 
k2 = 40 + A2+1 

} 

DielectricRegion //gridA-3 

{ 

QuseFlag = 0 
jl = 40 - A2 
kl = 40 + Al+1 
j2 = 40 - A2 
k2 = 40 + Al+1 

} 

DielectricRegion //gridA-8 

{ 

QuseFlag = 0 
jl = 40 + A1 
kl = 40 - A2+1 
|2 = 40 + A1 
k2 = 40 - A2+1 

} 

DielectricRegion //gridA-7 

{ 

QuseFlag = 0 
jl = 40 + A2 
kl = 40 - Al+1 
jjZ = 4 0 + A2 
k2 = 40 - Al+1 

} 

DielectricRegion //gridA-5 

{ 

QuseFlag = 0 
jl = 40 - A1 
kl = 40 - A2+1 
j2 = 40 - A1 
k2 = 40 - A2+1 

} 

DielectricRegion //gridA-6 

{ 

QuseFlag = 0 
jl = 40 - A2 
kl = 40 - Al+1 
j;2 = 40 - A2 
k2 = 40 - Al+1 

} 

//Regions for Anode grids 

DielectricRegion //anode grid 1 

{ 

QuseFlag = 0 
jl = 40 + FI 
kl = 40 + F2+1 
*§| = 40 + FI 
k2 = 40 + F2+1 
Secondary 
{ 

secondary = 1 
secSpecies = ions 
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iSpecies = ions 


I 

} 

DielectricRegion //anode grid 2 

{ 

QuseFlag = 0 
jl = 40 + F2 

kl = 40 + Fl+1 

j2 = 40 + F2 

k2 = 40 + Fl+1 

Secondary 

{ 


secondary = 1 
secSpecies = ions 
iSpecies = ions 


DielectricRegion //anode grid 4 

{ 

QuseFlag = 0 
jl = 40 - FI 
kl = 40 + F2+1 
Up = 40 - FI 
k2 = 40 + F2+1 
Secondary 
{ 

secondary = 1 
secSpecies = ions 
iSpecies = ions 

} 

} 

DielectricRegion //anode grid 3 

{ 

QuseFlag = 0 
jl = 40 - F2 
kl = 40 + Fl+1 
j;2 = 40 - F2 
k2 = 40 + Fl+1 
Secondary 
{ 


secondary = 1 
secSpecies = ions 
iSpecies = ions 


DielectricRegion //anode grid 8 

{ 


QuseFlag = 0 
jl = 40 + FI 
kl = 40 - F2+1 
jZ = 40 + FI 
k2 = 40 - F2+1 
Secondary 
{ 

secondary = 1 
secSpecies = ions 
iSpecies = ions 
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} 

DielectricRegion //anode grid 7 

{ 

QuseFlag = 0 
jl = 40 + F2 
kl = 40 - Fl+1 
j;2 = 40 + F2 
k2 = 40 - Fl+1 
Secondary 
{ 

secondary = '■% 
secSpecies = ions 
iSpecies = ions 

} 


DielectricRegion //anode grid 5 

{ 

QuseFlag = 0 
jl = 40 - FI 
kl = 40 - F2+1 
J2 = 40 - FI 
k2 = 40 - F2+1 
■Secondary 
{ 


secondary = 1 
secSpecies = ions 

iSpecies = ions 


DielectricRegion //anode grid 6 

{ 

QuseFlag = 0 
jl = 40 - F2 
kl = 40 - Fl+1 
j2 = 40 - F2 
k2 = 40 - Fl+1 
Secondary 
{ 


secondary = 1 
secSpecies = ions 
iSpecies = ions 


//add-in diagonal anode grids 


//Regions for Emitters deleted 


/////////////////////////////////////////End dieletric regions: Start 
Diagnostics: ////// 

Diagnostic 

{ 

Hi st Max = 1028 
Comb = 2 
Ave 0 
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ji = C 

kl = 40 

j2 = 80 

k2 = 40 

VarNarae = phi 

title = phi at equator 

xl_Label = equator position 

x2_T,abel; time 

} 

Diagnostic 

{ 

HistMax = 1028 
Comb = 2 
Ave = 0 
jl = 40 
kl = 40 
j2 = 40 
k2 = 40 
VarName = phi 

title = phi at center position 
xl_Label = time 
x2_Label = phi at center 

} 

Diagnostic 

{ 

HistMax = 1028 

Comb = 2 

Ave = 5 

jl - 0 

kl = 40 

j2 = 80 

k2 = 40 

VarName = rho 

title = rho at equator; 

xl_Label = equator position 

x2_Label = time 

} 

Diagnostic 

{ 

HiStMax = 1028 

Comb = 2 

Ave = 5 

jl = 30 

kl = 30 

j2 = 50 

k2 = 50 

VarName = rho 

title = rho in middle' 

xl_Label = x 

x2_Label = y 

} 

Diagnostic 

{ 

HistMax = 1028 
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Comb = .Jf : 

Ave = 5 
jl = 0 
kl = 45 
j2 = 80 
k2 = 45 
VarNarae = E2 

tiSJjfe = transverse acgelration along midlib® 

xl_Label = x 

x2_Label = E-transverse 

} 

} 
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